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Abstract 

Systems of differential-algebraic equations (DAEs) are generated routinely by simulation and mod¬ 
eling environments such as Modelica and MapleSim. Before a simulation starts and a numer¬ 
ical solution method is applied, some kind of structural analysis is performed to determine the 
structure and the index of a DAE. Structural analysis methods serve as a necessary preprocess¬ 
ing stage, and among them, Pantelides’s algorithm is widely used. Recently Pryce’s E-method is 
becoming increasingly popular, owing to its straightforward approach and capability of analyzing 
high-order systems. Both methods are equivalent in the sense that when one succeeds, producing a 
nonsingular system Jacobian, the other also succeeds, and the two give the same structural index. 

Although provably successful on fairly many problems of interest, the structural analysis meth¬ 
ods can fail on some simple, solvable DAEs and give incorrect structural information including the 
index. In this report, we focus on the E-method. We investigate its failures, and develop two 
symbolic-numeric conversion methods for converting a DAE, on which the E-method fails, to an 
equivalent problem on which this method succeeds. Aimed at making structural analysis meth¬ 
ods more reliable, our conversion methods exploit structural information of a DAE, and require a 
symbolic tool for their implementation. 
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Chapter 1 

Introduction 


We are interested in solving initial value problems in DAEs of the general form 

fi{t, the Xj and derivatives of them) =0, / = 1 : n, (1.1) 

where the Xj{t) are n state variables, and t is the time variable. The formulation (1.1) ineludes 
high-order systems and systems that are jointly nonlinear in leading derivatives. Moreover, (1.1) 
ineludes ordinary differential equations (ODEs) and purely algebraie systems. 

An important eharaeteristie of a DAE is its index. Generally, the index measures the diffieulty 
of solving a DAE numerieally. If a DAE is of index-1, then a general index-1 solver ean be used, 
e.g., DASSL [3], IDA of SUNDIALS [14], and MATLAB’s odelSs and ode2 3t. If a DAE is of 
high index, that is, index > 2, then we need a high-index DAE solver, e.g., RADAUS for DAEs of 
index < 3 [13] or DAETS for DAEs of any index [22]. We ean also use index reduetion teehniques 
to eonvert the original DAE to an index-1 problem [17, 19, 33], and then apply an index-1 solver. 

Structural analysis (SA) methods serve as a preproeessing stage to help determine the index. 
Among them is the Pantelides’s method [25], whieh is a graph-based algorithm that finds how many 
times eaeh equation needs to be differentiated. Pryee’s struetural analysis—the Signature method 
or L-method —is essentially equivalent to that of Pantelides [27], and in partieular eomputes the 
same structural index when both methods sueeeed. However, Pantelides’s algorithm ean only 
handle first-order systems, while Pryee’s ean be applied to (1.1) of any order and is generally 
easier to apply. 

This SA determines the struetural index, whieh is often the same as the differentiation index, 
the number of degrees of freedom, the variables and derivatives that need to be initialized, and the 
eonstraints of the DAE. We give the definition of the differentiation index in §2 and that of the 
struetural index in §3. 

Nedialkov and Pryee [20,21,22] use the E-method to analyze a DAE of the form (1.1), and 
solve it numerieally using Taylor series. On eaeh integration step, Taylor eoeffieients (TCs) for 
the solution are eomputed up to some order. These eoeffieients are eomputed in a stage-wise 
manner. This stage by stage solution seheme, also derived from the SA, indieates at eaeh stage 
whieh equations need to be solved and for whieh variables [24]. In [2, 12, 15], the E-method is 
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also applied to perform struetural analysis, and the resulting offset veetors are used to preseribe the 
eomputation of TCs. 

Although the E-method provably gives eorreet struetural information (ineluding index) on 
many DAEs of praetieal interest [27], it can fail—whence also Pantelides’s algorithm and other 
SA methods [34, 35] can fail—to find a DAE’s true structure, producing an identically singular 
system Jacobian. (See §3 for the definition of system Jacobian.) 

Scholz et al. [33] show that several simulation environments such as Dymola, OpenMod- 
ELICA and SimulationX all fail on a simple, solvable 4x4 linear constant coefficient DAE; we 
discuss this DAE in Example 4.18. Other examples where SA fails are the Campbell-Griepentrog 
Robot Arm [5] and the Ring Modulator [18]. When SA fails, the structural index usually underes¬ 
timates the differentiation index. In other cases, when SA produces a nonsingular system Jacobian, 
the structural index may overestimate the differentiation index [31]. We review in Appendix B how 
these DAEs in the early literature are handled so that SA reports the correct index. 

SA can fail if there are hidden symbolic cancellations in a DAE; this is the simplest case among 
SA’s failures. However, SA can fail in a more obscure way. In this case, it is difficult to understand 
the causes of such failures and to provide fixes to the formulation of the problem. Such deficiencies 
can pose limitations to the application of SA, as it becomes unreliable. Our goal is to construct 
methods that convert automatically a system on which SA fails into an equivalent form on which 
it succeeds. This report is devoted to developing such methods. 

It is organized as follows. Chapter 2 overviews work that has been done to date. Chapter 3 
summarizes the E-method and gives definitions and tools that are needed for our theoretical devel¬ 
opment. The problem of SA’s failures on some DAEs is described in Chapter 4. In Chapters 5 and 
6, we develop two methods, the linear combination method and the expression substitution method, 
respectively. We show in Chapter 7 how to apply our methods on several examples. Chapter 8 gives 
conclusions and indicates several research directions. 


Chapter 2 

Background 


The index of a DAE is an important eoneept in DAE theory. There are various definitions of 
an index: differentiation index [4, 9, 10], geometrie index [30, 32], struetural index [7, 25, 27], 
perturbation index [13], tractability index [ 11 ], and strangeness index [16]. 

The most commonly used index is the differentiation index', we refer to it as d-index or Vd- The 
following definition is from [1, p. 236]. 

Definition 2.1 Consider a general form of a first-order DAE 

F(t,x,x')=0, (2.1) 

where dE/dx' may be singular. The differentiation index along a solution x{t) is the minimum 
number of differentiations of the system that would be required to solve x' uniquely in terms of x 
and t, that is, to define an ODE for x. Thus this index is defined in terms of the overdetermined 
system 


dF 

dt 


F(t,x,x') =0, 
(t,x,x',x") =0, 


dPF 


t,x,x,--- ,x 


(p+1) 


0 


to be the smallest integer p so that x! in (2.1) can be solved for in terms ofx and t. 


If a DAE (1.1) is of high-order, then one can introduce additional variables to reduce the order 
of the system so that it is still in the general form ( 2 . 1 ). 

We give a definition for solution of a DAE. 

Definition 2.2 An n-vector valued function x(t), defined on a time interval I C M, is a solution of 
(1.1), if (t ,x(t)) satisfies fi = 0, i= 1 : n, pointwise for all t E I: that is, eve/ 7 /;■ vanishes on I. 
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ReiBig et al. [31] claim that a DAE of d-index 1 may have arbitrarily high structural index. 
They construct a class of linear constant coefficient DAEs in some specific form. On these DAEs 
of d-index 1, Pantelides’s algorithm performs a high number of iterations and differentiations, and 
obtains a high structural index that far exceeds the d-index 1. A simple 3x3 linear electrical circuit 
example is also presented: choosing a specific node as the ground node results in a DAE of d-index 
1 , but of structural index 2 . 

Pryce [27] shows that, if the E-method succeeds, then the structural index V 5 is always an 
upper bound on the d-index. This implies that, if the structural index computed by the E-method 
is smaller than the d-index, then the method must fail; otherwise we would have a statement that 
contradicts to the above Definition 2.1. Pryce also shows that the E-method succeeds on one of 
ReiBig’s DAEs and produces a nonsingular system Jacobian [27]. His method also produces the 
same high structural index as does Pantelides’s. 

In [26], Pryce shows that the E-method fails on the index-5 Campbell-Griepentrog Robot Arm 
DAE—the SA produces an identically singular Jacobian. He then provides a remedy: identify 
the common subexpressions in the problem, introduce extra variables, and substitute them for 
those subexpressions. The resulting equivalent problem is an enlarged one, where the E-method 
succeeds and reports the correct structural index 5. Pryce introduces the term structure-revealing 
to conjecture that a nonsingular system Jacobian might be an effect of DAE formulation, but not 
of DAE’s inherent nature. 

Choudhry et al. [ 6 ] propose a method called symbolic numeric index analysis (SNIA). Their 
method can accurately detect symbolic cancellation of variables that appear linearly in equations, 
and therefore can deal with linear constant coefficient systems. Eor general nonlinear DAEs, SNIA 
provides a correct result in some cases, but not all. Eurthermore, it is limited to order-1 systems, and 
it cannot handle complex expression substitution and symbolic cancellations, such as (vcosy)' — 
y cosy. Eor the general case, their method does not derive from the original problem an equivalent 
one that has the correct index. 

Scholz et al. [33] are interested in a class of DAEs called coupled systems. In their case, a 
coupled system is composed by coupling two semi-explicit d-index 1 systems. They show that the 
E-method succeeds if and only if the coupled system is again of d-index 1. As a consequence, if 
the coupled system is of high index, SA methods must fail. They develop a structural-algebraic 
approach to deal with such coupled systems. They differentiate a linear combination of certain al¬ 
gebraic equations that contribute to singularity, append the resulting equations, and replace certain 
derivatives with newly introduced variables. They use this regularization process to convert the 
regular coupled system to a d-index 1 problem, on which SA succeeds with nonsingular Jacobian. 


Chapter 3 


Summary of Pryce’s structural 
analysis 


We call this SA [27] the E-method, because it constructs for (1.1) an n x n signature matrix E = 
{Oij) such that 



the order of the highest order derivative to which xj occurs in /,■; or 
—oo if Xj does not occur in ft. 


(3.1) 


A transversal E is a set of n positions (/, j) with one entry in each row and each column. The 
sum of entries Oij over T, or Y.{ij)eT is called the value T, written Val(r). We seek a highest- 
value transversal (HVT) that gives this sum the largest value. We call this number the value of the 
signature matrix, written Val(E). 

We give a definition for a DAE’s structural posed-ness. 

Definition 3.1 We say that a DAE is structurally well-posed (SWP) if its Val{'L) is finite. That is, 
all entries in a HVT are finite, or equivalently, there exists some finite transversal. Otherwise, if 
Valiff) — then we say a DAE is structurally ill-posed (SIP). 

For a SWP DAE, we find equation and variable ojfsets c and d, respectively, which are non¬ 
negative integer n-vectors satisfying 

Ci > 0; dj — Ci > Oij for all i, j with equality on a HVT. (3.2) 

An equality dj — ci = Oij on some HVT also holds on all HVTs [29]. We refer to c and d satisfying 
(3.2) as valid offsets. They are not unique, but there exists unique c and d that are the smallest 
component-wise valid offsets. We refer to them as canonical offsets. 

The structural index is defined by 



max, Ci -1-1 if = 0 for some j, or 
max,c, otherwise. 
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Critical to the success of this method is the nonsingularity of the DAE’s nxn system Jacobian 
matrix J = (J,;), where 

T dfi f if J —c,-= and 

Jo = ^ 

dxj ^ [ 0 otherwise. 

Note that J = J(c,d) depends on the choice of valid offsets c,d, which satisfy (3.2). That is, 
using different valid offsets, one may obtain different system Jacobians. However, they all have the 
same determinant; see Theorem 4.15. For all the examples in this report, we shall use canonical 
offsets and the system Jacobian derived from them. 

We can use E and c, d to determine a solution scheme for computing derivatives of the solution 
to (1.1). They are computed in stages 


+ 1 ,..., 0,1,... where kd = — max 

j 


At each stage we solve equations 

0 = for all i such that c, + A: > 0 (3.4) 

for derivatives 

for all j such that dj + k>0 (3.5) 

using the previously found 


for all j such that 0 <r < dj-\-k. 

We refer to [24] for more details on this solution scheme; see also Example 3.2. 

Throughout this report, for brevity, we write “derivatives of xf instead of “vj and derivatives 
of it”—derivatives of a variable v include v itself as the case 1 = 0. 

If the solution scheme (3. 4-3. 5) can be carried out up to stage k = 0, and the derivatives of 
each variable xj can be uniquely determined up to order dj, then we say the solution scheme and 
the SA succeed. The system Jacobian is nonsingular at a point 


(f Y, ro 

It, , A25 . . . 5 ^ 


-Ui), 



(3.6) 


and there exists a unique solution through this point [20, 27, 29]. We say the DAE is locally 
solvable, and call (3.6) a consistent point, if derivatives xp^ do not occur jointly linearly in fp\ 
In the linear case, a consistent point is 


r. ro 

t , A1, . . . , Aj 5 A25 • 


,,^2 


(^ 2 - 1 ). 


^ri’) • 


jAn 1 , 


(3.7) 


For a more rigorous discussion of a consistent point, we refer the readers to [20, 24, 29]. 
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To perform a numerical check for SA’s success, or a success check for short, we attempt to 
compute numerically a consistent point at which J is nonsingular up to roundoff: we provide an 
appropriate set of derivatives of x/s and follow the solution scheme (3.4-3.5) for stages k = kj:0. 
This set of derivatives is the set of initial values for a DAE initial value problem, and a minimal 
set of derivatives required for initial values is discussed in [29]. 

When SA succeeds, the structural index is an upper bound for the differentiation index, and 
often they are the same: < V 5 [27]. Also, the number of degrees of freedom (DOF) is 

DOF = Val(E) = ^J,-^c,= ^ a.y. 

./ '■ ('j)er 

We say the solution scheme and S A fails, if we cannot determine uniquely a consistent point 
using the solution scheme defined by (3.4-3.5) —otherwise said, we cannot follow the solution 
scheme up to stage k = 0 and find a consistent point at which J is nonsingular. In our experience, 
in the failure case usually Vd > Vs, but not always, and the true number of DOF is overestimated 
by Val(E). This is discussed in Examples 4.7, 4.9, 4.18, 4.19. 

We illustrate the above concepts using the following example. 

Example 3.2 The simple pendulum DAE (Pend) in Cartesian coordinates is 

0 = /i = x" -\-xX 

^ = f 2 =y" +y^-g (3.8) 

0 = f3=x^W-L\ 

Here the state variables are x, y, A; g is gravity, and L > 0 is the length of the pendulum. 

The signature matrix and system Jacobian of this DAE are 



X 


A 

Ci 


X 

3^ 

A 

/l 

'2* 

— 

o' 

0 

/l 

' 1 

0 

X 

E= /2 

— 

2 

0* 

0 

and J = -^2 

0 

1 

y 

h 

. 0 

0* 

— 

2 

h 

. 

2y 

0 . 


dj 2 2 0 


We write E in a signature tableau: a HVT is marked by •; — denotes — 0 °; the canonical offsets 
c, d are annotated on the right of E and at the bottom of it, respectively. 

The structural index is 

Vs = maxc,- + 1 = C 3 + 1 = 3, 

i 

which is the same as the d-index. The number of degrees of freedom is 

DOF = ^J,-^c, = 2. 

.7 ' 
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stage k solve for using previously found 

-2 0 = /3 Ji;,y 

-1 0 = /' /,y 

>0 0 = y*+2)^3;(fc+2)^;|^(fe) ^{<k+2)^y{<k+2)^y^{<k) 


Table 3.1: Solution seheme for (3.8) 


Sinee the derivatives j = 1,2,3, that is, x'',y",X, oeeur jointly linearly in (3.8), a eonsis- 
tent point is given by ,y,y'). If we evaluate J at this point, then 

det(J) = -2{x^+y^) = -2L^ ^ 0 

(beeause +y^ = l} by fi = 0) and SA succeeds [27]. The solution scheme is in Table 3.1. The 
notation is short for z,z',... 

For brevity, in the following chapters, when we give a system of equations, we write down 

• the signature matrix, 

• a HVT in it (marked by •), 

• the canonical offsets c, d, 

• positions (/, ]) where dj — Ci > Oij > 0 (marked by ^), and 

• the accompanying system Jacobian. 

When we present a SA result, we omit the words 

“the signature matrix and system Jacobian are in the following.” 

Provided there is a finite HVT in E, we also show the value of the signature matrix and the deter¬ 
minant of the system Jacobian—Val(E) and det(J). For instance, after giving (3.8), we simply put 
E with Val(E) attached, and J with det(J) at the bottom. 



X 

J 

A 

Ci 


X 

3^ 

A 

/l 

'2* 

— 

o' 

0 


' 1 

0 

X 

E= /2 

— 

2 

0* 

0 

J= ^2 

0 

1 

y 

h 

. 0 

0* 

— 

2 

h 

. 2.^ 

2y 

0 . 


det(J) = -2L2 


2 


2 


0 


Val(E) = 2 
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Similarly, if we write the signature matrix of a system as E, then we write correspondingly 
the canonical offsets as c, d, and the Jacobian as J. Throughout this report, we shall show DAE 
problems for which our conversion methods are suitable. These methods critically depend on the 
SA results. 



Chapter 4 

Structural analysis’s failure 


In this chapter, we investigate how SA fails on some DAEs. That is, SA produces a singular system 
Jacobian, and the problem is solvable. In §4.1, we give definitions for (a) a structural zero in the 
system Jacobian, and (b) a structurally singular DAE, where the system Jacobian is identically 
singular. In §4.2 we identify two types of SA’s failure. 

4.1 Success check 

To perform a success check for SA on a SWP DAE, we attempt to evaluate the system Jacobian J 
in (3.3). If a point (3.6) satisfies the solution scheme (3.4-3.5) at stages k = + 1,..., 0, and J 

is nonsingular, then SA succeeds. 

In the definitions that follow, we let A be an n x n matrix function. 

Definition 4.1 An {i,j) position is a structural zero of A if Aij is identically 0; otherwise it is a 
structural nonzero. 

Definition 4.2 [20] Matrix A is structurally singular if every B G with Bij = 0 inA’s struc¬ 

tural zero positions, is singular — equivalently, if every transversal of A contains a structural zero. 
Otherwise A is structurally nonsingular. 

Definition 4.3 Matrix A is identically singular, if its determinant is identically 0,' otherwise it is 
generically nonsingular. 

Eor a matrix function, being structurally singular is a special case of being identically singular; 
see Example 4.4 below. 

Example 4.4 Consider the following three matrix functions of variables x and y: 


Ai = 


X 

0 




X X 


A y' 

A 2 = 


, and A 3 = 


y y. 


A 


12 
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Ai is identically singular because det(Ai) = 0. It is also structurally singular, since every 
B G with B 21 = B 22 = 0 is singular. Here, (2,1) and (2,2) are structural zero positions of A, 
and each transversal in A contains a structural zero. 

A 2 is also identically singular, as det(A 2 ) = xy — xy = 0. It is structurally nonsingular, since a 
transversal does not contain a structural zero. 

A 3 is structurally nonsingular. It is generically nonsingular, since det(A 3 ) = x^ —y^ is not 
identically zero. A 3 is singular only when x = ±y. 

In the following, we denote (1.1) by and define two concepts for it: 

• a structural zero in the system Jacobian J, and 

• a structurally singular DAE. 

Let J be the set of index-pairs 

jr={( 7,/)|7 = l:n,/eN}. (4.1) 

Given an n-vector function x = x(t) that is sufficiently smooth (but not necessarily a solution of 
J”), let 

= I 

For a finite subset J of J, we define a |7|-vector xy whose components are x^p as ( 7 , /) ranges over 
J. (The ordering of these components does not matter.) 

Now we denote a DAE as J^. We define the derivative set of as 

derset(J^) = { ( 7 , Z) | x^ occurs in (4.2) 

Then the derivatives occurring in can be denoted concisely as X(jerset( 7 ^)- 

By a value point we mean a G M x that contains values for t and values for the 

derivative symbols in . 

Example 4.5 In the simple pendulum DAE (3.8), the state variables x,y, A are xi,X 2 ,X 3 . Eet L = 5 
and g = 9.8. Then 

derset(.F) = { (1,0), (1,2), (2,0), (2,2), (3,0) }. 

A possible value point can be 


*3 (L -^1 At A 2 A 2 As) (A 53 , 3,4,1.6,1), 


which satisfies /i and fs but not / 2 . 
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Similarly, we define the derivative set o/J: 
derset(J) = { (y, 1) \ oeeurs in J}. 

From (3.3), a derivative oeeurring in J must also oeeur in T', but not viee versa. For example, 
in Pend, xf',y",X do not appear in J, and derset(J) = {(1,0), (2,0) }; ef. Example 3.2. The 
derivative set of J is a subset of that of iF: derset(J) C derset(J^). 

Definition 4.6 An (i, j) position is a struetural zero o/J, ifiij is identically zero at all value points 
G M X that satisfy 0 or more equations from 

0 = fp\ m>0, i=l:n. (4.3) 

Otherwise, (z, j) is a struetural nonzero. 

For the present purpose, we do not require the DAE to have a unique solution, or even any 
solution. That is, we do not eonsider existenee and uniqueness of the DAE at this stage, while 
identifying struetural zeros of J and the singularity of J diseussed below. 

Reeall (3.3) that defines J. If dj — Ci > Oij, then J,y = 0 and thus position (z,y) is a struetural 
zero in J. The eonverse is not true; see Example 4.7. 

Example 4.7 Consider an artificially modified simple pendulum DAE. We multiply the first equa¬ 
tion /i by + y^ — if and obtain 


0 = fi = {x"+xX){x^+y^-L^) 

0 = f 2 = y"+y^ -g (4.4) 

0 = h=x^+y^-L\ 



X 

y 

A 

C/ 


X 

J 

A 

/l 

'2* 

0 

0 

0 

/i 

P 

0 

xjl 

E= /2 

— 

2 

0 * 

0 

J= ^2 

0 

1 

y 

h 

. 0 

0 * 

— 

2 

h 

. 2-^ 

2y 

0 . 

dj 

2 

2 

0 

Val(E) = 2 

det(J) = 

-2ll{x^+y 


In J, /i =x^ +y‘^ —L^. To decide which entries in J are structural zeros, we notice the following. 

• If we evaluate J at some random then ji is not identically equal zero. Hence positions 
(/i,x) and (/i. A) are not identical zeros. 
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• If we evaluate J at some ^ that satisfies 

At =/3 = 0, 

then aeeording to Definition 4.6, positions (/i,x) and (/i, A) are structural zeros of J. 

We give a definition for structural regularity of a DAE. 

Definition 4.8 A DAE is structurally singular if J is identically singular at all value points ^ G 
M X that satisfy 0 or more equations from (4.3). Otherwise the DAE is structurally 

nonsingular, or structurally regular. 

Example 4.9 In the previous example, positions (/i,x) and (/i, A) are structural zeros of J at any 
point that satisfies /s = 0. By Definition 4.8, (4.4) is structurally singular. 

In fact, it can be shown that a solution of Pend is a solution to (4.4), but not vice versa. 

Example 4.10 Consider the DAE in [1, p. 235, Example 9.2], written in (1.1) form: 

0 = /i = -y'l +y3 

0 = /2 = y2(i-y2) (4.5) 

0 = /3 = JD2+y3(l -y2) 



yi 

^2 

3^3 

Ci 



3^2 y3 

/l 

' r 

— 

o' 

0 

h 

'-1 

0 1 

E= ■/2 

— 

0 * 

— 

0 

J= /2 

0 

l-2y2 0 

h 

.0 

0 

0 *. 

0 

h 

. 0 

yi-y3 1-3^2. 

dj 

1 

0 

0 

Val(E) = 1 

det(J) = 

-(1 -2y2)(l -y2) 


SA gives V 5 = 1, and det(J) depends solely ony 2 - Erom /2 = 0, either y 2 = 0 ory 2 = 1- To examine 
ifj is nonsingular, we consider each of the following two cases. 

• If y 2 = 0, then det(J) = — 1 and SA succeeds. In this case (4.5) is of d-index 1. 

• If y 2 = 1, then det(J) = 0 and SA fails. This failure comes as no surprise because (4.5) is 
now of d-index 2 and SA underestimates its index; see the discussion in §2. 

Remark 4.11 Eor a structurally ill-posed (SIP) DAE, there does not exist a finite transversal in its 
E—every transversal in E contains at least one — 0 °. In this case, there exists no valid offsets c, d, 
not to mention a system Jacobian that depends on these offsets. In contrast, a structurally singular 
DAE has valid offsets and a system Jacobian that is identically singular. Hereby we distinguish the 
difference between a SIP DAE and a structurally singular DAE. 
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Suppose J is generically nonsingular. If J is singular when evaluated at a point along a solution, 
then we say the DAE is locally unsolvable at this point, and we call it a singularity point. See 
Example 4.12. 


Example 4.12 [8] Consider 

0 = /i = -x' + y 
0 = /2 = a: + cos(t)y. 


(4.6) 



X 

3^ 

Ci 


X 


ft 

'l* 

o' 

0 

h 

'-1 

1 

h 

_0 

0 *. 

0 


. 0 

cos(t)_ 

dj 

1 

0 

Val(E) = 1 

det(J) = 

— cos(t) 


Since det(J) is generically nonzero, (4.6) is structurally nonsingular. We can integrate this 
problem from t = 0 with any consistent initial value (v(0),y(0)) = (vo,yo)> and the problem is 
index-1 (both differentiation and structural indices) as long as det(J) ^ 0. However, J is singular 
att = tk= {k+l/2)n,k = 0,l,.... Hence, we say the DAE has a singularity point at tk- 


4.2 Identifying structural analysis’s failure 

We give below a definition for the true highest-order derivative (HOD) of a variable xj in a function 
u. 

Definition 4.13 The true HOD of xj in u is 

( 1 highest order derivatives ofxj on which u truly depends; or 

I ^ —oo ifu does not depend on any derivative of Xj (including xj). 

By “truly” we mean that, if r = a [xj, u) > then u is not a constant with respect to x^p. Eor 
example, u = x'-\-cos^x”-\- sin^x" = x'-|- 1 truly depends on x' but not x”, resulting in a (x,u) = 1. 

In practice, however, we usually find the formal HOD of xj in u, denoted by a (xy, u) , instead of 
the true HOD. By “formal” we mean the dependence of an expression (or function) on a derivative 
without symbolic simplifications. Eor example, u = xf cos^x" -|- sin^x" formally depends on x" 
and hence a (x, u) = 2, while w = x' -I- 1 and a{x,u) = 1 . 

We denote also Oij = d {xj.fi) corresponding to Ojj. The daets and daesa codes implement 
[21, Algorithm 4.1 (Signature matrix)] for finding formal dij. 

Since the formal dependence is also used in [21, §4], we can adopt the rules in [21, Eemma 
4.1], which indicate how to propagate the formal HOD in an expression. The most useful rules are: 
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• if a variable v is a purely algebraie funetion of a set U of variables u, then 

a (xi,v) = max a (xj,u), 

' ueU ' 


(4.7) 


and 

• if V = where p > 0, then 

a {xj,v) = a {xj,u) + p. (4.8) 

These rules are proved in [21], to whieh we refer for more details. We illustrate the rules in 
Example 4.14. 

Example 4.14 Let u = {X 1 X 2 )' — x\x 2 . Applying (4.7) and (4.8), we derive the formal HOD of xi 
in u: 


o{xi,u) = max{ a (xi, (X 1 X 2 )'), o [xi,x[x 2 ) } 

= max{ o{xi,x\X 2 ) + 1 , max{ o {xi,x'i) ,o{x\,X 2 ) } } 

= max{max{ o (xi,xi), o (xi,X 2 ) } + l,max { 1,-°°} } 

= max| max| 0 , — 0 ° } + 1 , 1 } 

= max| 0 + 1 , 1 } 

= 1 . 

Similarly o {x 2 ,u) = 1. Simplifying u = {xiX 2 )' —x\x 2 results in m = xix^. Henee, the true HOD 
of x\ in M is (7 (xi,m) = 0, and that of X 2 in w is a (x 2 ,m) = 1. 

When sueh a hidden symbolic cancellation oeeurs, o [xj, u) ean overestimate the true o [xj, u). 
If u is an equation //, then the formal HOD o {xj,fi) may not be the true Oij. We write Ojj = 
o (xj-ifi) eorresponding to Oij = o (xj,fi). We call the matrix E = (dij) the “formal” signature 
matrix. Also, let c, d be any valid offsets for E, and let J be the resulting Jacobian defined by (3.3) 
with E and c, d. 

If Oij > Oij, then // does not depend truly on x^p^\ That is, fi is a constant with respect to x^p^\ 
Then = 0, and {i,j) is a structural zero in J. Due to such cancellations, J has more structural 
zeros than J does, and hence J is more likely to be structurally singular. It is also possible that the 
DAE itself is structurally ill posed. 

Since Oij for all i, j = I : n, we can write E > E meaning “elementwise greater or equal”. 

We define the essential sparsity pattern Sess of E to be the union of the HVTs of E. That 
is, the set of all {i,j) positions that lie on any HVT. We give two theorems below, which are 
Theorems 5.1 and 5.2 in [20]. In the following, we use the term “offset vector” to refer to the 
vector (c,d) = (ci,... ,c„, Ji,... ,<i„). 
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Theorem 4.15 Suppose that a valid offset vector (c, d) for E gives a nonsingular J as defined by 
(3.3) at some consistent point. Then every valid offset vector gives a nonsingular J (not necessarily 
the same as J) at this point. All resulting J, including J, are equal on Sess, tind all have the same 
determinant det(J) = det(J). 

By “equal on Sess” we mean Jij = J,y for all (/, j) e S^ss- 

Theorem 4.16 Assume that J, resulting from E and a valid offset vector (c, d), is genetically 
nonsingular. Let (c, d) be a valid offset vector for the formal signature matrix E, and let J be the 
Jacobian resulting from E and (c, d). In exact arithmetic, one of the following two alternatives 
must occur: 

(i) VafiL) = Val{L). Then every HVT ofL is a HVT ofL, andc, d are valid offsets for E. Conse¬ 
quently, J is also generically nonsingular. 

(ii) Val(’L) > VaKY). Then J is structurally singular. 

Theorem 4.16 shows that J, resulting from E > E and a valid offset veetor (c, d), is either 

(1) nonsingular, and SA is using valid, but not neeessarily eanonieal, offsets for the true E; or 

(2) strueturally singular, and SA fails due to symbolie eaneellations, in a way that may be deteeted. 

In the latter ease, this failure may be avoided by performing symbolie simplifieation on some 
or all of the /;’s. However, “no elever symbolie manipulation ean overeome the hidden eaneella- 
tion problem, beeause the task of determining whether some expression is exaetly zero is known 
to be undeeidable in any algebra elosed under the basie arithmetie operations together with the 
exponential funetion” [20]. 

Example 4.17 Consider 

D = f\ = (xy)'-x'y-xy +2x + y-3 
0 = f 2 =xLy- 2 . 



X 

J 

Ci 


X 

J 

ft 

' r 

1 ’ 

0 

ft 

' 0 

o' 


. 0 

0*. 

1 

A 

1 

1 


dj 1 1 Val(E) = l det(J)=0 


Here, the signature matrix and Jaeobian are the formal ones. Sinee det(J) = 0, SA fails. Simplify¬ 
ing /i to /i =2x-\-y — 3 reveals that (4.9) is a simple linear algebraie system: 

0 = /i = 2x + y-3 
0 = f2^xLy-2. 






CHAPTER 4. STRUCTURAL ANALYSIS’S FAILURE 


19 



X 


Ci 


X 


h 

'o* 

o' 

0 

/l 

' 2 

r 


. 0 

0*. 

0 ■>= 

h 

1 

1 

dj 

0 

0 

Val(E) = 0 

det(J) = 

1 


Another kind of SA’s failure oeeurs when J is not strueturally singular, but is identieally singu¬ 
lar. Examples 4.18 and 4.19 illustrate this ease. 

Example 4.18 Consider the coupled DAE from [33] ^ : 


0 = /i = -x[+X 3 +bi{t) 

0 = /2 = -X2+X4+b2{t) 

0 = / 3 = X2+X3+X4+Ci{t) 

0 = /4 = -Xi +X3 +X4+C2{t). 



Xi 

X 2 

■t:3 

X4 

Ci 


Xl 

•tt2 

-t:3 

X4 

/i 

'r 

— 

0 

— 

0 

/i 

'-1 

0 

1 

o' 

/2 

— 

1* 

— 

0 

0 

fi 

0 

-1 

0 

1 

h 

— 

0 

0* 

0 

0 

/3 

0 

0 

1 

1 

h 

_ 0 

— 

0 

0*. 

0 

/4 

. 0 

0 

1 

1_ 


dj \ 1 0 0 Val(E) = 2 det(J) = 0 


This DAE is of d-index 3, while SA finds structural index 1 and singular J. Hence SA fails. 
Example 4.19 In the following DAE, SA reports the correct d-index 2 but still fails. 

0 = /i = -x\-X3+XiLX2+gl{t) 

0 = /2 = -X 2 -X 3 +Xi -hX2 -hX3 +X4+g2{t) 

0 = /3 =X2+X3+g3{t) 

0 = f4=Xi-X4+g4{t) 


'We consider this DAE with parameters b — e — 1, ai — a 2 — 5 — 1, and 7 = —1. In [33] superscripts are used 
as indices, while we use subscripts instead. We also change the (original) equation names gugi to / 3 ,/ 4 , and the 
(original) variable names y \,72 to X 3 ,X 4 . 
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Xl 

JC 2 

■^3 

X4 

Cf 


Xl 

JC 2 


X4 

/l 

’l* 

0 

1 

— 

0 

fl 

'-1 

0 

-1 

o' 

fl 


1 * 

1 

0 

0 

fl 

0 

-1 

-1 

1 

/3 

— 

0 

0 * 

— 

1 

h 

0 

-1 

-1 

0 

/4 

[0 

— 

— 

0 *. 

0 

h 

. 0 

0 

0 

1 _ 


dj 1 1 1 0 Val(E) = 2 det(J) = 0 


Using the solution scheme derived from the SA result, we would try to solve at stage k = 0 the 
linear system 0 = /i,/ 2 ,/ 3,/4 for where the matrix is J. Since it is singular, the 

solution scheme fails in solving (4.11) at this stage; see Table 4.1. 


stage k solve for using comment 

— 1 0 = /3 X 2 ,X'i — initialize xi 

0 0 = /i,/ 2 ,/ 3,/4 x^,X2,V3,4:4 x\,X 2 ,X 2 , J is singular; solution Scheme fails 


Table 4.1: Solution scheme for (4.11) 

Now we replace /2 by f 2 — f 2 + to obtain 

0 = /i = -x\-X2,+XiLX2Lgl{t) 

0 = /2 +/3 = /2 = +X2 +X3 +X4 +g2(0 +.?3(0 

0 = f3=X2+X3+g2{t) 



0 

II 

= Xl - 

X 4 +g 4 (t). 







Xl 

JC 2 

X 3 

X4 

Ci 


Xl 

JC 2 

X 3 

X4 

fl 

' 1 

0 

1 * 

— 

0 

fl 

'-1 

0 

-1 

0 

72 

0 * 

0 

0 

0 

1 

fl 

1 

1 

1 

1 

/3 

— 

0 * 

0 

— 

1 

/3 

0 

1 

1 

0 

/4 

. 0 

— 

— 

0 *. 

1 

/4 

_ 1 

0 

0 

-1 


dj \ \ \ \ Val(E) = 1 det(J) = 2 


(4.12) 


The solution scheme succeeds; see Table 4.2. The resulting DAE (4.12) is of structural index 
V 5 = 1 , which equals the differentiation index. 
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At stage k = 0, we solve 0 = for x'-^,x! 2 ,x! 2 ,x'^ using x\,X 2 ,X'i,x^. Sinee /2 = 

fi /s ’ need f” to find these first-order derivatives. Therefore, the original DAE (4.11) is of 
differentiation index 2. 

Note that by setting fi = f 2 ~ /s we ean immediately reeover the original system. It can be 
easily verified that a vector function 

X{t) = (xi(t),X2(t),X3(t),X4(t))^ 

that satisfies (4.12) also satisfies (4.11), and vice versa. We explain in §5 how this conversion 
makes SA succeed. 


stage k solve for using comment 

-1 0 = 72,/3,/4 Xi,X2,X3,X4 - 

0 0 = /i,72:/3:/4 jci,X 2 ,X 3 ,X 4 J is nonsingular; 

solution scheme succeeds 

Table 4.2: Solution scheme for (4.12) 


In Examples 4.18 and 4.19, J is not structurally singular, but is still identically singular. No 
symbolic cancellation occurs in the equations therein. Therefore, this kind of failure is more diffi¬ 
cult to detect and remedy, and we wish to find techniques to deal with such failures. 

We call our techniques conversion methods, and describe them in the upcoming chapters. We 
wish to convert a structurally singular DAE into a structurally nonsingular problem, provided some 
conditions are satisfied and allow us to perform a conversion step. The original DAE and the 
converted one are equivalent in the sense that they have (at least locally) the same solution set. We 
shall also elaborate on this equivalence issue. 



Chapter 5 

The linear combination method 


In this chapter we introduce the linear combination method, or the LC method for short. We present 
in §5.1 some preliminary lemmas. Then we describe in §5.2 how to perform a conversion step. In 
§5.3 we give definitions and results about equivalence of DAEs and address how equivalence is 
related to the LC method. 

For simplicity, throughout this report, we consider only the second type of SA’s failures de¬ 
scribed in §4.2: “singular” means identically singular but not structurally singular. Based on this 
assumption, symbolic cancellations are not considered an issue that makes the E-method fail. 


5.1 Preliminary lemmas 

Lemma 5.1 (Griewank’s Lemma) [21, Lemma 5.1] Let v be a function oft, Xj’s and derivatives 
of them (j = \ : n). Denote v^p^ = d^v/dt^, where p > 0. If o {x q, then 

dv _ dv' 
dxf ~ dx^^^^^' 

Hence 

dv dv' dv^P^ 

dxf ~ dx^^^^'^ ~ dx^^^P^' 

Lemma 5.2 Let E and Lbe nxn signature matrices. Assume VaHT) is finite, c, d are valid offsets 
for E, and Oij < dj — cifor all z, j—l:n. If a HVT in E contains a position (z, j) where Oij < dj — cu 
then Val{L) < Val(L). 

Proof. Let T be a HVT in E. Then 

Val(E) = ^ o,j <j^dj-j^a = Val(E). □ (5.2) 

(f;)er j=i z=i 
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Corollary 5.3 For a row index I, let 

{ Oij = Oij for all i 7 ^ I and all j, and 

Oij<dj-ci for all]. 

Then Val{E) < Val{E). 

Proof Since Gij < dj — ci for all j, the intersection of a HVT in E with positions in row I is a 
position (/, r) with Oir < — c/. By Lemma 5.2, Val(E) < Val(E). □ 

This lemma shows that, if we replaee a row / in E with a row with entries less than dj — ci for 
eaeh column j, then the value of this signature matrix deereases. 


5.2 Conversion step 

Given a SWP DAE of the form (1.1), assume that we apply the E-method and obtain a singular 
system Jacobian J. We seek a reformulation of this DAE so that the system Jaeobian J of the new 
DAE may be generieally nonsingular. We denote by E and E the signature matrices of the original 
DAE and this new DAE, respectively. Denote by c, d the valid offsets for E. 

We deseribe below how to perform a eonversion step using a linear eombination (EC) of equa¬ 
tions. We call this conversion technique the LC conversion method, or simply the LC method. The 
main result from this conversion is that, under eertain eonditions, we can obtain an equation in a 
row, say I, sueh that xj occurs in this row of order < dj — ci for all j. Henee by Corollary 5.3, 
Val(E) < Val(E). 

We assume n>2. Eet w be a nonzero vector function from the null space of J^. Here, J and u 
are eonsidered as funetions of t, xj’s and appropriate derivatives of them. 

Denote by I{u) the set of indices for whieh the ith component of u is not identieally zero 

I{u) = {i\my^0}, (5.3) 

and let 

6{u) = min c/. (5.4) 

i(zl{u) 

Sinee u is nonzero and J is identieally singular, I{u) has at least two elements. Otherwise J has a 
row of identieal zeros and is structurally singular. 

Remark 5.4 We eonsider u in its simplest form in the sense that its elements do not have a eom- 
mon factor comprising t, xfs, or/and derivatives of them. Eor instance, in Example 4.18, we do 
not use u = {Q,Q,x!-^,—x\Y though = 0, but use u = (0,0,1, — 1)^. 

Also, we do not eonsider u with any fraetions. Eor example, we use u = (0, 0,x'^,x\X2Y instead 
of ( 0 , 0 ,Jl; 7 ^X 2 (Vl)^l)^. 


CHAPTER 5. THE LINEAR COMBINATION METHOD 


24 


The sufficient condition for applying the LC method is the following: for a nonzero u G ker(J^), 


o (xj, u) < dj — 6{u) for all j = \ :n 


(5.5) 


If this condition is satisfied, then we can perform a conversion step. We explain this “sufficiency” 
in Remark 5.8. 

Denote by L{u) C I(u) the set of indices I such that the /th component of u is not an identical 
zero and c/ = 0(m) = min;g/(„) cf. 

L{u) = [l e I{u) I c/ = 0 (m) }. (5.6) 


From (5.4), there exists at least one I G I{u) such that c; = 0{u), so L{u) ^ 0. 
We choose an / e L{u) and replace fi by 

7, = E 

iCiI(u) 


(5.7) 


We refer to (5.7) as a conversion step using the LC method and to the resulting DAE as a converted 
DAE. Critical for the success of the EC method is the following lemma. 

Theorem 5.5 For a SWP DAE with identically singular J, let u be a nonzero n-vector such that 
fu = 0. If 

o {xj,u) < dj — G{u) for all j = I \ n 


and we replace fi by fi in (5.7), then the converted DAE has E with VaKT) < VaKY). 

Eirst, we illustrate with an example how to perform a conversion step, and then we prove this 
lemma. Since u is fixed during a conversion step, for brevity we write I{u), 9{u), and L{u) as 7, 9, 
and L, respectively ^ 


Example 5.6 Consider 

0 = /l = -x[ +X 3 
0 = f2 = —Xj + Xa 

. (5.8) 

0 = /3 =F{xi,X 2 ) 

0 = /4 = XsFx^ [xi ,X 2 )+ X4Fx^ [xi , X 2 ) + G{xi , X 2 ). 

Here [xi ,X 2 ) = (9F(xi,X 2 ) / dxi , and similarly we write {xi , X 2 ), Gx^ [xi , X 2 ), and Gx^ {x\ ,X 2 ). 


^This set L is not to be confused with the constant L in the pendulum-related DAEs. 
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Xl 

JC2 

■^3 

X 4 

C, 


Xl 

71:2 

•^3 

X 4 

/l 

’l* 

— 

0 

— 

0 

fl 

'-1 
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fl 

— 

1 
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0* 

0 

fl 

0 

-1 

0 

1 

h 

0 

0* 

— 

— 

1 

h 



0 

0 

/4 

[0 

0 

0* 

0. 

0 

h 

. 0 

0 

Fx, 

Fx2 \ 

dj 

1 

1 

0 

0 

Val(E) = 1 


det(J) 

= 0 



Because of singular J, the SA fails. It reports structural index 2, but the differentiation index is 3. 
If we take u = (F^j, I, —l)^, then = 0. 

We illustrate (5. 3-5. 6): 

/={/| m,^0} = {1,2,3,4}, 

9 = mine, = 0, 
iel 

L={lel\ci = e=0} = {l,2,4}. 

Then we check if condition (5.5) holds: 

o{xi,u)= 0 <1 = ^1 — 0 , 
o{x 2 ,u)= 0 < 1=^2 — 0 , 

(7 (x 3 ,m) = —oo < 0 = J 3 — 0 , and 
O (x 4 ,m) = —00 < 0 = J 4 — 0. 

Hence o [xj, u) < dj — 0 for all j. 

Using (5.7) gives 

iel iel 

= Fxifl+ Fx2f2 + f 2 ,— f4 

= Fxj ( -x\ + X3) + Fc2 (-4 +M) + {F{xi,X2))' - (x 3 Fxj + X4FX2 + G) 

= -XjFt, +X3Ft, -x'2Fx2 Fx^Fx^ Fx\Fx^ +X 2 FX 2 -x^Fx^ -x^^Fx^ - G 
= -G. 

Now, with 0=0, 

o{xiJ) = 0 < 1 = Ji-0, 

(y{x2,f) = 0 < 1 = J 2 - 0 , 

o [x 2 ,,f) ——00 <0 = d^ —6, and 
a {x4,f) = —00 < 0 = J 4 — 0 . 
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That is, <7 < dj —6 for all j. 

For each I e L = {1,2,4}, assuming ui ^ 0, we can replace // by /; = /. We show in the 
following the three possible converted DAEs, each with Val(E) = 0 and generically nonsingular J. 


• 1 = 1 : 


0 = 7l = -G(jCi,JC2 ) 

0 = /2 = -4 +X4 
0 = /3 =F(xi,X2) 

0 = /4 = (xi ,X 2 )+ X4F_,2 (xi , X 2 ) + G(xi , X 2 ) 
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— 

— 
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Fx2 

0 
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/4 
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0 

0* 

0. 

0 

f4 

0 

0 

Fx, 

Fx2. 

'^2 

1 

1 

0 

0 

Val(E) = 0 

da{i) = Fx^{Fx,Gx^ 

-Fx2 

Gx,) 


When Ml = Fii 7 0 and F^, determinant is nonzero and the SA succeeds. 


• 1 = 2 : 


0 = /l = -x[ +JC3 

0 = 72 = -G( 3 Ci, .3:2) 

0 = /3 =F{xi,X 2 ) 

0 = /4 = v:3F;tj (.ri, 3 C 2 ) + X 4 FX 2 (.ri, X 2 ) + G{xi , X 2 ) 


(5.10) 
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0 
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0 
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1 

1 

0 

0 

Val(E) = 0 

= Fx2{Fx,Gx2 

-Fx2 

Gx,) 


Similarly, the SA succeeds when M 2 = 7 0 and Fx^ Gx^ 7 • 
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• / = 4: 


0 = /i = -x[ +^3 
0 = /2 = -X2 +X4 
0 = /3 =E(xi,X 2 ) 

0 = f4^-G(xi,X2) 
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1 
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0 

Val(E) = 0 

det(J) = - 

-Fxi Gx 2 

+ Fx2 

Gxi 


In this case, SA’s success requires only -fxi Gx2 FxjGxx ■ 


Using the LC method, we obtain three converted DAEs from (5.8): (5.9), (5.10), and (5.11). 
However, it is not guaranteed that all converted DAEs and the original DAE have exactly the same 
solution sets. We will cover this equivalence issue in §5.3. 

Now we prove Eemma 5.5. 

Proof. We show first that 


Oij = o {xjJi) < dj - Cl for all j = I : n. 


Since 0 = min,g/c,, c,- — 0 > 0 for i G I. By (3.2), a {xj^fi) = (Jij < dj — ci. Applying 
Griewank’s Eemma (5.1) to (3.3), with q = Ci — 0, gives 


Sfi SfC”'' 




for i G I. 


(5.12) 
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Then for all 7 = 1 : n, 


i=l 


= 7 Ui 


dfi 




E"' 


SllisiUif 


(o-e; 


dxf 


dfi 


dx] 


(dj-e) 



using = 0 

iel 





using (5.12) 


using 0 (xy, m) 


using (5.7). 


(5.13) 


(dj^e) 


This shows that fi does not truly depend on Xj , that is, 
Gij = o < dj — 0 = dj — Cl for all j = I : n. 

By Corollary 5.3, Val(E) < Val(E). 


□ 


Remark 5.7 We name this method “LC” beeause of the following eonsiderations. The veetor u 

T id—6) 

from the null spaee of y that satisfies (5.5) does not eomprise the leading derivatives x - ^ for 
all j in equation f^‘^‘ for all i G /. We consider here each ui as a “constant” in 

7, = E“.#''”. 

iel 

and fi as a “linear combination” of equations 

Remark 5.8 If o [xj, u) =dj — 0 for some j, then Val(E) is not guaranteed < Val(E). In this case, 
we cannot swap the sum and the differentiation operator in (5.13). Therefore, we cannot prove 

dfi/dx^^^ = 0. Then, in E, it can happen that 

'Oij = o (xj-ifi) = dj —9 = dj — Cl for some 7 , 
giving < instead of strictly < in (5.2). 

However, if o (xj, u) < dj — 9 holds for 7 from a particular set, we can still achieve Val(E) < 
Val(E). We leave this investigation for future work, consider the condition (5.5) sufficient for now, 
and require it to be satisfied for the LC method. 

If M is a constant vector, then o [xj, u) = — 0 ° for every Xj, and the condition (5.5) is automat¬ 
ically satisfied. In this case we do not need to check it. We illustrate this in the next example. 
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Example 5.9 Consider 

0 = /i =Xi+tX2+t^X3+gi{t) 

0 = f 2 =x\+ + g 2 [t) 

0 = /s = x" -\-tX2-\- + ^ 3(0 



Xl 

2 C2 

X3 

Ci 


JCl a:2 

JC3 

/l 

’o* 

0 

o' 

2 

/i 

’ 1 t 

r2' 

E= /2 

1 

r 

1 

1 

J= ^2 

1 t 

r2 

h 

_ 2 

2 

2 *. 

0 

/3 

. 1 t 

2^1 

dj 

2 

2 

2 

Val(E) = 

3 

det(J) 

= 0 

= (-1,1, 

0)^ 

,J^u 

= 0. 

Using (5.3 

-5.6) gives 



/={1,2 


e = 

C2 = 

1 , and 

L={2}. 




Since w is a constant vector, condition (5.5) is satisfied. We replace /2 by 
= +U 2 f 2 ^^'’ 

= -/l'+/2 

= — (xi + tX2 + t^X3 — .?l) + {x'l + tx'2 + t'^x'3 +.^2) 

= -v:2-2r2C3-/i+^2- 

The converted DAE is 

0 = /i =Xi+tX 2 Et^X 3 +gi 
0 = f 2 = -X 2 - 2 tX 3 -g\ +g 2 
0 = /3 = x'( + tX 2 + 21 ^X 3 + g 3 . 



Xl 

2^2 

2C3 

Ci 


2Cl 

3^2 

X3 

/i 

'o* 

0 

o' 

2 

/i 

' 1 

t 


E= /2 

— 

0 * 

0 

2 J = 

72 

0 

-1 

-2t 

h 

_ 2 

2 

2*. 

0 

/3 

_ 1 

t 

2t\ 

dj 

2 

2 

2 

Val(E) = 2 


det(J) = 

= -r2 


lit > i/e for a suitable £ depending on the machine precision, then J is computably nonsingular. 
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5.3 Equivalent DAEs 


First, we give a definition for equivalent DAEs. 

Definition 5.10 Let and denote two DAEs. They are equivalent (on some interval for t) if a 


solution of(F is a solution to and vice versa. 

In the following context, we denote by T the original DAE with equations //, i = \ : n, and 
singular Jacobian J. After a conversion step using the EC method, we obtain a (converted) DAE, 
denoted by IF, with equations /,-, z = 1 : n, and Jacobian J, which may be singular. 

Theorem 5.11 After a conversion step using the LC method, DAEs IF and IF are equivalent on 
some real time interval I/or t, ifui Ofor all t G I. 

Proof. Eet a solution of IF, over some interval I C M, be a vector-valued function 

X(t) = {xi{t),...,Xn{t)Y 

that satisfies (1.1) for all t G I. 

We denote the vector used in the EC method by zz = (zzi,..., zz„), where its zth component is of 
the form 



If u is defined at (/x(t)) for all t G I, then 


fi = 'Lu,fY and /. = /.forzV/ 


vanish at (/x(t)), and thus x(t) is a solution to IF. 

Conversely, assume that x(t) is a solution of IF on I. If u is defined at (/x(t)) for all t G I and 
ui then 

// = -(//- E and /■=7,forzV/ (5-14) 

vanish at (/x(t)), and thus x{t) is a solution to IF. 

By Definition 5.10, and IF are equivalent. □ 

Remark 5.12 We can see from (5.14) that, if we have a choice for I, it is desirable to choose it 
such that ui is identically nonzero, e.g., a nonzero constant, jcj -|- 1, or 2 -|-cos^X 3 . In this case, IF 
and IF are always equivalent—we do not need to check the equivalence condition zz/ 7 0 when we 
solve IF. 


Example 5.13 In Example 5.6, case 1=1 [resp. 1 = 2] requires Fx, 7 0 [resp. F;c 2 7 0] to recover 
the original DAE (5.8) from (5.9) [resp. from (5.10)]. However, for case / = 4, zz 4 = 1 is a nonzero 
constant for any t. Therefore this choice is more desirable than the other two. 
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Below we define an ill-posed DAE using the struetural posedness defined in the DAESA papers 
[23,29]. 

Definition 5.14 A DAE is ill posed if it has an equivalent DAE that is structurally ill-posed (SIP); 
otherwise it is well posed. 

Example 5.15 Consider problem (4.4). Using 0 = fs = -hy^ — L^, we reduce /i to the trivial 
0 = /i = 0. This is just performing a simple substitution, and is not applying the LC method. The 
signature matrix 


E = 


7i 

/2 

/3 


X y X 

-20 
0 0 - 


(5.15) 


does not have a finite HVT, so the resulting DAE is SIP. Hence, by Definition 5.14 , the original 
SWP DAE (4.4) is ill posed. 

Corollary 5.16 If a structurally well-posed DAE can be converted, by the LC method, to an equiv¬ 
alent DAE that is structurally ill-posed, then the original DAE is ill posed. 

Proof. This follows from Theorem 5.11 and Definition 5.14. □ 

Example 5.17 Consider the following SWP DAE 

0 = /i = y”+yX -\-yX' 

^ = f 2 = y"+yX-g (5.16) 

0 = h=x^W-L\ 



X 

y 

X 

Ci 


X 


A 

fl 

— 

3 

1 *' 

0 

fl 

'o 

1 

y 

E= /2 

— 

2 * 

0 

1 

J= 72 

0 

1 


h 

0 * 

0 

— 

0 

h 

_2x 

0 

0 . 

dj 

0 

3 

1 

Val(E) = 

3 

det(J) 

= 0 

Eor M = (1, 

-1,0)^, fu 

= 0. Using 

(5. 3-5. 6) gives 



'={1.2}. 

e = 

Cl = 

= 0, and 

Z,= {1|, 
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Since m is a constant vector, condition (5.5) is satisfied. We replace /i by 

7l = /l - /2 = (/' +3^^') - (/ +3^^ -^)' = 0. 

The signature matrix of the resulting problem is exactly (5.15). Hence, by Corollary 5.16, (5.16) 
is ill posed. 

If the Jacobian of the converted DAE is still singular, we may be able to apply the LC method 
iteratively, provided condition (5.5) is satisfied on each iteration. Since after each conversion step 
we reduce the value of the signature matrix by at least 1 , the number of iterations does not exceed 
Val(E), where E is for the original DAE. We use Example 5.18 to show how we can iterate with 
the EC method. 

Example 5.18 We construct the following (artificial) ModPendA DAE from Pend (3.8): 

0 = A = f2,-\-f \ + (x + xX ) 

0 = B = fi +A" = x” +xX + {x^ +y^-L^ + (x" +xA)')" (5.17) 

0 = C = f 2 +A'" = y" +yX -g + [x^+y^ -L^ + {x" ExX)')'”. 



X 

3^ 

A 

Ci 


X y 

A 

A 

’3* 

0 

1 ' 

3 

A 

’l 2y 

X 

e0 = ^ 

5 

2 * 

3 

1 

j0_ 5 

1 2y 

X 

C 

.6 

3 

4*. 

0 

C 

. 1 2y 

x_ 

dj 

6 

3 

4 

Val(EO) 

= 9 

det(jO) 

= 0 


Here, a superscript denotes an iteration number, not a power. We show how to recover the simple 
pendulum problem. 

We find a vector in ker((J®)^): w® = (—1,1,0)^. Then 

7® ={1,2}, 00 = 1, and L^ = { 2 ]. 

We replace the second equation B by 

_a(3-i) +5 = _A" + (A" + /i) = /i = x" + xA. 

The converted DAE is 

0 = A =x^+y^-L^ + {x +xX)' 

0 = /i = x^^ +xA 

0 = C = y' +yX -g + [x^ +y^ -L} + [x' +xX)')'”. 
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X 

y 

A 

c, 


X y 

A 

A 

'3* 

0 

r 

3 

A 

1 2 y 

X 

E^ = /i 

2 

— 

0 * 

4 

Ji= /i 

1 0 

X 

C 

.6 

3* 

4_ 

0 

C 

. 1 2 y 

X. 

dj 

6 

3 

4 

Val(Ei) 

= 6 

det(J^) 

= 0 


Although Val(E^) = 6 < 9 = Val(E®), matrix is still singular. If = (—1,0,1)^, then 
= 0. This gives 

/^ = {1,3}, 0^=0, and l1 = {3}. 

We replaee the third equation C by 

_a(3”0) + c = -A'" + (/2 + A'") = /2 = / + yA - g. 

The eonverted DAE is 

0 = A =x^+y^-L^ + ix' +xXy 
0 = /i = x!' +xX 
0 = /2 =/+yA -g. 



X 

y 

A 

Ci 


X 

y 

A 

A 

’3* 

0 

r 

0 

A 

’l 

0 

X 

I? = fi 

2 

— 

0 * 


/i 

1 

0 

X 

h 

— 

2 * 

0 . 

0 

/2 

.0 

1 

0 . 

dj 

3 

2 

1 

Val(E2) = 5 


det(J^) 

= C 


We have Val(E^) = 5 < 6 = Val(E^), but is still singular. We find u = (1, — 1,0)^ sueh that 

(J 2 )rj^ 2_0 

/2 = { 1 , 2 }, 02 ^ 0 , and L^ = {\]. 

Replaeing the first equation A by 

A -/( = {h+f[)-f[ =h=x^ W-L\ 

we reeover /i,/ 2,/3 from (5.17). This is exaetly the DAE Pend (3.8), with Val(E) = 2 and 
det(J) = — 2L^; ef. Example 3.2. 

Sinee eaeh u in every eonversion iteration is a eonstant veetor, eaeh ui we piek is a nonzero 
eonstant. By Remark 5.12, the original DAE (5.17) and Pend are always equivalent. Henee, we 
ean solve (5.17) by simply solving Pend. 













Chapter 6 

The expression substitution method 


We develop in this ehapter the expression substitution eonversion method. In §6.1, we introduee 
some notation. We deseribe in §6.2 how to perform a eonversion step using this method and address 
in §6.3 the equivalenee issue. 


6.1 Preliminaries 

A eonversion using the LC method seeks a row in E, replaees the eorresponding equation by a 
linear eombination of existing equations, and eonstruets a new DAE with a signature matrix of a 
smaller value. Inspired by the LC method, our goal is to develop a eonversion method that seeks 
a eolumn in E, performs a ehange of eertain variables, and eonstruets a new DAE with E sueh that 
Val(E) ^^al^E^. 3A^e refer to this approaeh as the ajcpt'assxof'i suixstitutioyi tnathod^ 

or the ES method. 

Again, we start from a SWP DAE with a signature matrix E, offsets c, d, and identieally singu¬ 
lar Jaeobian J. To start our analysis, we give some notation below. 

Let M be a veetor funetion from the null spaee of J, that is, Ju = 0. Denote by L{u) the set of 
indiees j for whieh the jth eomponent of u is not identieally zero 

L{u) = { j \ Uj ^0}, (6.1) 

and denote s{u) by the number of elements in L(m): 

s{u) = \L{u)\. (6.2) 

Note that s >2. Otherwise J has a eolumn that is identieally the zero veetor, and henee J is 
strueturally singular. 

Let 

I{u) = { / I — Ci = Oij for some j G L{u) }. (6.3) 


34 
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Denote also 


C(M) = maxc,-. (6.4) 

ieliu) 

Now we illustrate (6.1-6.4). 

Example 6.1 Consider 

0 = fl=Xl+ -^ 2^2 + ^1 (0 

/ 9 .X (6.5) 

0 = f 2 =Xi +X 2 X 2 +X2+g2it)- 



Xl 

■^2 

Ci 


Xl 

JC2 

/l 

'l* 

2 ' 

0 

/i 

-p 

-11X2 

A 

. 0 

1* 

1 


1 

X2 _ 


dj 1 2 Val(E) = 2 det(J) = 0 

In J, /i = Using u = (jC 2 , —1)^ for which Ju = 0, ( 6 .1-6.4) become 

L(m) = { 1 , 2 }, 5(m) = |L(m)| = 2 , 

Am) = {1,2}, (6.6) 

and C(u) = max c, = C 2 = 1. 

iEl(u) 

We show later how the ES method works on this problem. 

T 

Remark 6.2 Assume that we apply the LC method to (6.5). First, we find u = (1,/i) from 
ker( J^). Using the notation in the LC method, we find / ={l,2},0=O, and k=\. Since 

C7(xi,m) = C7(xi,jU) = 1 = 1 - 0 = - 0, 

the condition (5.5) is not satisfied. After a conversion step, the resulting DAE is still structurally 
singular with Val(E) = Val(E) = 2 and identically singular J. See §7.3 for more details. 


6.2 A conversion step using expression substitution 

We can perform a conversion step using the ES method, if the following conditions hold for some 
nonzero u such that Jw = 0: 

i \ / f dj-C{u)-\ ify gL(m) 

<3 [Xj, U) < < 

1 dj — C{u) otherwise 


(6.7) 
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and 


dj—C{u) > 0 for all j G L{u) 


( 6 . 8 ) 


We call (6.7) and (6.8) the sufficient conditions for applying the ES method. 

Picking an / G L, we show below how to perform a conversion step. Since we use the same u 
throughout the following analysis, we omit the argument u and simply write L, s, I, and C. 
Without loss of generality, assume that the nonzero entries of u are in its first s positions: 

u = 

Then L= { j \ uj = {1,..., 5 }, where s = \L\. 

We introduce ^ variables y 1 ,..., and let 


yj=^j 


[dj-o 


^jJd,-C) 

A/ 

ui 




yi = x'l 


(di-C) 


for; gL\{Z}, 
for 7 = Z. 


(6.9) 


(The condition (6.8) guarantees that the order of xj, j G L, in (6.9) is nonnegative.) 
Written in matrix form, (6.9) is 


y\ ■ 


■ 1 

-Ul/ui 



r 1 

Xi 

yi 

= 


i 



Jdi-C) 

ys . 



-Us/Ul 

1 _ 


(dl~C) 


This s X s square matrix is nonsingular with determinant 1. 
We write the first part of (6.9) as 




(dj-C) 


= yj+-N 

Ul 


jJd,-C) 


for jeL\{l}. 


By (6.4), Ci < C for all Z G I. Differentiating (6.10) C — c,- > 0 times yields 


'^K-oyc-c,) 


yj+—N 

Ul 


^jJd,-C) 


(C-Ci) 


In each f with Z G /, we replace every with Oij = dj — ci and j G L \ { Z } by 




( 6 . 10 ) 


Denote by /, the equations resulting from these substitutions. For Z ^ 7, we set fi = f. 
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From (6.9), we introduce the equations that prescribe the substitutions: 



for; gL\{/} 


for j = 1. 


( 6 . 11 ) 


We append these equations to the /,’s and construct an augmented DAE that comprises 

equations Jj, gi,... in variables xi,... .Xn, yu.. .ys- 

We write the signature matrix and system Jacobian of this converted problem as E and J, respec¬ 
tively. 

Example 6.3 For (6.5), assume we pick I = 2. Since 


L\{/} = {1.2}\{2} = {1}. 


we introduce y; = yi. Since 

= 1, ^2 = 2, Cl = 0, C = C 2 = 1, and w = (jC 2 , —1)^, 


(6.10) becomes 



In /i, we replace x"f^ = Vj by 


(yi -X 2 X 2 )(‘^ = (yi -X 2 X 2 )(^ = (yi -X2X2)' = y'j -^2 -X2X2. 


Similarly, in /2 we replace x\ by yi —^ 2 X 2 . Taking these substitutions into account and appending 
gi and g 2 , we obtain 


0 = = XI +^;-(3'1-^24)'-X24' 

= xi+c->'i+^?+gi(t) 


0 = /2 = ( j^l --^ 24 ) +-^24 +^2 +,? 2(0 

= yi+V2+g2(0 


0 = gi = -yi-hxi-hX 2 X 2 
0 = g 2 ^ -yi+M- 


( 6 . 12 ) 
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/l 

X\ 

■ 0 

JC2 

1 * 

1 

3^2 

Ci 

0 

/i 

Xl 

' 1 

JC2 

Ix^OL 

3^1 

- a 

3^2 

0 

h 

— 

0 

0 * 

— 

1 

fi 

0 

2X2 

1 

0 

g\ 

0 * 

1 

0 

— 

J = 

0 

gi 

1 

X2 

0 

0 

g2 

_ — 

0 

— 

0 *. 

0 

g2 

. 0 

0 

0 

-1 

dj 

0 

1 

1 

0 

Val(E) = 1 

det(J) 

= Jf2- 

2a{x2 +^2 

) 


Here a = e . If det(J) ^ 0, then SA sueeeeds on (6.12). 

Our aim is to show that Val(E) < Val(E) after a eonversion step, provided that the suffieient 
eonditions (6.7) and (6.8) hold. Before proving this inequality, we state two lemmas related to the 
strueture of E. 


Lemma 6.4 Let c = (ci,..., c„) and d = (Ji,..., dn) be the valid offsets for E. Let c and d be the 
two {n + s)-vectors defined as 



if] = \:n 

if j = : n + 5, and 

ifi = \ '.n 

ifi = n-\- \ '.n-\-s. 


Then E is of the form in Figure 6.1. 


(6.13) 

(6.14) 


The proof is rather teehnieal, and we present it in Appendix A. 

Lemma 6.5 After a conversion step using the ES method, ValfL) < ValfL). 

Proof. If E does not have a finite HVT, then Val(E) = —oo, while the original DAE is SWP with a 
finite Val(E). Hence Val(E) < Val(E) holds. 
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Xl ■■■ X/_i 

Xl 

Xl+\ ' ' ' Xg 

-^j+i 

Xn 

y\ ■■■ yi-i 

yi 

yi+i ■■■ ys 

/l 


< 


< 


< 

— CO 

< 

fn 







— CO 


gl 

< 

= 

< 

< 


0 


— oo 

gl 


= 


— CO • • • 

— CO 


0 



< 


< . 

< 


— oo 



gs 


= 

= 





0 

d] 

di • ■ ■ <i/_i 

di 

d-i+i ■■■ ds 

ds+i 

dn 

C ■■■ C 

c 

c ■■■ c 


Figure 6.1: The form of E for the eonverted DAE using the ES method. The <, <, = mean 
the relations between Otj and dj — cu respeetively. Eor example, an [ij) position with < has 
~ D- 

On the other hand, if E has a finite HVT T so that Val(E) > —oo, then 
Val(i)= £ a,j 

sinee dj — Ci > Oij for all i,j= 1 : n + s 

using (6.13) and (6.14), 


(y)er 

n+s n+s 

= 

j=l i=l 


'^dj + sC j - { '^Ci + sC 
0=1 


0=1 


£j,-£c, = Val(E). 
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We prove in the following that Val(E) = Val(E) leads to a eontradietion. Assume this equality 
holds. Then there exists a transversal T of (n + s) positions in E sueh that 

dj — Ci = Oij > —oo for all {i,j) G T. (6.15) 

The eolumn eorresponding to yi has only one finite entry „+/ = 0, and therefore 
(n + Z,n + /) G T. Consider (Zi, 1),.. ., G T. Since (n + Z,n + Z) G T, row numbers Zi,...,Zi 
take values among 

1,2, ...,n, n + l,...,n + Z — l,n + Z+l,...,n + j'. (6.16) 

In (6.16) only ^ — 1 numbers are greater than n. Hence at least one of these row numbers is 
among 1 : n. That is, there exists (Z,., r) G T with 1 < Z^ < n and 1 < r < 5. This entry is in 
E(1 : n, 1 : 5 )^ in Figure 6.1, so dr — Ci, < which contradicts to our assumption in (6.15). 
Therefore Val(E) < Val(E). □ 


Remark 6.6 We give several remarks about the ES method. 

• After a conversion step, we perform symbolic simplifications on f^, for Z G /. By doing this 

we ensure that the ^'^’s for 7 G L = {1,. .. , 5 } disappear from these equations. That is, 
o < dj — Ci for 7 = 1 : 5 and Z G I. 

• Clearly, yi appears only in g/. We mark down the positions in T on E, and then remove row 
n + Z (corresponding to gi) and column n + Z (corresponding to y/). Because (n + Z, n + Z) G T, 
the remaining marked positions still form a HVT T in the resulting (n + 5 — 1) x (n + 5 — 1) 
signature matrix E. 

Since On+i,n+l = Val(E) = Val(E). The purpose to use gi and yi in the above proof and 
analysis is for our convenience. In practice, we can exclude gi and yi in the resulting DAE. 
For consistency, after removing g/ and y;, we still use E and T to denote the signature matrix 
and system Jacobian, respectively, for the resulting DAE. See Example 6.7. 

• If some derivativefor Z = 1 : n and 7 G L \ | Z }, appears implicitly in an expression 

in fi, then we need to write this expression into a form in which x^^'' appears explicitly. 
See Example 6 . 8 . 


'Using MATLAB notation. 
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Example 6.7 Removing g 2 and y 2 in (6.12) results in a DAE with the following signature matrix 
and Jacobian 



Xl 

2^2 

yi 

Ci 


Xl 

2^2 

yi 

/l 

■ 0 

1 

1 *' 

0 

/l 

' 1 

2^2 ct 

-a 

E= /2 

— 

0 * 

0 

1 

J= fl 

0 

2X2 

1 

^1 

.0* 

1 

o_ 

0 

gl 

_ 1 

2C2 

0 


dj 0 1 1 Val(E) = 1 det(J) = —X 2 + 2a(x2+X 2 ) 


Here a = 

Example 6.8 Suppose that x'l' appears implicitly in (sin2xj)" in /i, and that the ES method finds 
L={l,2}, / = 2, ^ 1=^2 = 3, and C = ci = 0. 


We want to replace vf ‘ = x'/' by 


3^1 + 


Ml (d2-C)\^^ 


U2 


U] 

yi H- 

U2 


fff 

2 • 


To make x'[' appear explicitly in /i, we expand (sin2xj)" and write 2x'['co^x\ — A{x'[Ysmx\ in¬ 
stead. Now we can perform the substitution forx'". 


6.3 Equivalence for the ES method 

We discuss here the equivalence for the ES method. Our approach below is similar to that for 
deriving the equivalence for the EC method. 

We denote by T the original DAE with equations fu i = \ : n, and a singular Jacobian J. 
After a conversion step using the ES method, we obtain a (converted) DAE T with equations 
i=l\n-\-s, and a Jacobian J, which may be singular. Here = gj, j = l ■. s. 

Assume that 

T 

x(t) = (xi(t),...,x„(t)) 

is a solution of T on some real time interval I C M. That is, every //, i=\-.n, vanishes at (t, x(t)). 
Assume also u is defined at (t,x(t)) for all t G I. We can substitute x(t) in (6.9) to find 

y(0 = 

such that every = gj, 7 = 1 : 5 , in (6.11) vanishes at (t,x(t),y(t)). Using (6.9), we perform 
substitutions in /,, i G I, to obtain /,. We let /,• = fi for i ^ 7. Since these substitutions do not 
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change the value of eaeh equation, eaeh /,■ also vanishes at (t,x(t),y(t)). Therefore (x(t),y(t)) is 
a solution to 

Conversely, assume that (x(t),y(t)) is a solution of on I C M. Assume also that u is defined 
at (t, x(t), y{t )) for all t G I. Note here u depends merely on t and x(t). Sinee ui is a denominator in 
each gj in (6.11), this solution requires ui{t) 7^ 0 on I. Given that eaeh gj vanishes at (t,x(t),y(t)), 
from (6.9) we have 





V 



^ / 




( 9 ) 


(?) 


jeL\{Z} 

j = i. 


where q >0. Substituting the expressions on the right-hand side for the derivatives of yj in each 
fi recovers ft and does not ehange its value. Therefore, eaeh fi also vanishes at (t,x(t),y(t)), or 
simply (t,x(t)) sinee y(t) does not appear in fi. Then x(t) is a solution to T. 

The above diseussion gives 

Lemma 6.9 After a conversion step using the ES method, DAEs T and T are equivalent ifui 7 ^ 0 
for all t G I. 

Again, if we have a choice for Z, it is desirable to choose one (whenever possible) such that ui 
is identically nonzero. In this case, T and T are always equivalent and we do not need to cheek 
U[f^0 when we solve J^. 


Example 6.10 In (6.5), assume we piek Z = I. Then (6.10) beeomes 
x '2 = =yiE = y2 -xilx2. 

U\ 

Here we use 

^Zi = I, d 2 = 2, C=I, and u = {x 2 ,—l)^. 

Then we 


substitute for in 
(y2--^iA2)' x '2 f\ 
yi-X\lx2 4 /2 

The equations gj derived from ( 6 . 1 1) are 

0 = gi = -yi +xi 
0 = g2 = -y2 +X'2+Xi/X2. 
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As / = 1, we can remove gi and yi, append equation g = g 2 , and obtain the resulting DAE 
0 = f^=xi+ 

= XI + e--^i-^2/2-4^! A 2 +x'i ^ 

0 = /2 = XI EX2{y2-X\lx2) +X 2 +g2{t) 

= X2y2+^2+S2{.t) 

Q = g= -y2+X2+Xi/x2. 



XI 

X2 

3^2 

Ci 



XI 


3^2 

/l 

0 

1 

1 *' 

0 

fi 

’1 

-X 2 / 3 /X 2 

-X 1 / 3 /X 2 

-X2P 

72 

— 

0 * 

0 

1 

J= f2 


0 

2X2+y2 

X2 

g 

. 0 * 

1 

0 . 

0 

g 


IA 2 

1 

0 

dj 

0 

1 

1 

Val(E) = 

1 det(J) = 

= -X2 + /3 (2x2 + y2 + ^2 

-X 1 /X 2 ) 


In J, /3 = exp(—X 2 y 2 “ 42 Xi/x 2 ). If det(J) y 0, then SA succeeds and gives structural index V 5 = 2. 
Here Val(E) = 1 < 2 = Val(E). 

However, the original DAE and the resulting one are equivalent only if mi = X 2 4 0 on some 
time interval I. In practice, it is more desirable to choose 1 = 2 since w/ = — 1 is identically nonzero; 
see also Example 6.3. 






Chapter 7 

Examples 


In this chapter, we illustrate how to apply the LC method and the ES method to several strueturally 
singular DAEs. When a eonversion method succeeds, we obtain an equivalent strueturally regular 
DAE with a nonsingular system Jaeobian. 

In §7.1, we apply both eonversion methods to the 4x4 linear eonstant eoeffieient (eoupled) 
DAE (4.10). The EC method sueeeeds in eonverting this problem to a strueturally regular DAE 
in two iterations, redueing the value of the the signature matrix by 2. In eontrast, the ES method 
reduees the value of the signature matrix by 1 in the first iteration. In the seeond iteration, the 
eondition for applying the ES method is not satisfied, and henee it eannot be applied further. 

In §7.2, we illustrate both methods on an artifieially eomplieated problem ModPendB derived 
from the simple pendulum DAE Pend (3.8). We show in §7.2.1 how the ES method sueeeeds in 
eonverting this problem to a strueturally regular DAE, whieh has a relatively simple strueture. In 
§7.2.2, the EC method is applied, but yields a eonsiderably more eomplieated result. 

In §7.3, we address Remark 6.2 in more detail: the eondition for applying the EC method is 
not satisfied for (6.5). If we perform a eonversion step, then the value of the signature matrix is not 
guaranteed to deerease. 


7.1 A simple coupled DAE 

Reeall the 4x4 linear eonstant eoeffieient (eoupled) DAE (4.10): 

'0 = /i = -x\+x^ + bi{t) 

^ 0 = /2 = -4 +-^4 + b2{t) 

0 = /3= X2+X3+X4+Ci(t) 

^ 0 = /4 = -Xi +X3 +X4 +C2(t). 


(7.1) 


44 
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/l 

Xl 

'r 

3C2 

-3^3 

0 

X 4 

Ci 

0 

fi 

X\ 

'-1 

3C2 

0 

^(^3 

1 

JC4 

o' 

/2 

— 

1 * 

— 

0 

0 

fi 

jO = 

/3 

0 

-1 

0 

1 

/3 

— 

0 

0 * 

0 

0 

0 

0 

1 

1 

/4 

. 0 

— 

0 

0 *. 

0 

/4 

. 0 

0 

1 

1 _ 

dj 

1 

1 

0 

0 

Val(EO) 

= 2 

det(jO) = 

0 



We use^ to denote the original problem. We let E® and denote the signature matrix and 
Jaeobian of the original problem, respeetively. 

7.1.1 LC method 

We show how to apply the LC method to this problem, and finally obtain an equivalent strueturally 
regular DAE on whieh SA sueeeeds. 

Let = (0,0, — 1,1)^. Then = 0. Using (5. 3-5. 6) gives 

/0={3,4}, 00 = 0, and l0 = {3,4}. 

We ehoose Z® = 3 G 7® and replaee fs by 

73 = “3/3 + W4/4 = -/3 +/4 = -Xl -X 2 -C 1 (t) + C2(0- 

The eonverted DAE is 

' 0 = fi = -x'l +X3 + bi (t) 

0 = f2 = -X2+X4 + b2(t) 

■ < - 

0 = fj = -XI -X2-Ci(t)+C2(t) 

, 0 = f4 = -Xi+X3+X4+C2(t). 




X] 

- 3:2 

• 3:3 

X 4 

Ci 


X\ X 2 

- 3:3 

X 4 

/l 


1* 

— 

0 

— 

0 

fl 

■ -1 0 

1 

0 ■ 

fl 


— 

1 

— 

0* 

0 

fl 

0 -1 

0 

1 

f3 


0 

0* 

— 

— 

1 

f3 

-1 -1 

0 

0 

/4 


0 

— 

0* 

0 . 

0 

h 

0 0 

1 

1. 



1 

1 

0 

0 

Val(Ei) 

= 1 

det(J^) 

= 0 



' We use a superscript to mean an iteration number, not a power. 
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Since is still singular, we apply again the LC method. Let = (—1, — 1,1,1)^. Then 
= 0. Now 

/I = { 1,2 ,3,4}, 0i=O, and l1 = {1,2,4}. 

We ehoose / ^ = 1 G and replaee /i by 

7i = “i/i + M 2 / 2 +“ 3 / 3 +“ 4/4 

= —fl —fl+Zs+fA 

= - [-x\ Ex-i+biit)] - [-x' 2 +x^ + b 2 {t)\ + [-Xi -X 2 -Ci{t) +C 2 {t)]' 

+ [—X\ + 3 C 3 -\-X 4 -\- C2(t)] 

= -JCi -blit) -b 2 {t) -c'lit) +4(0 +C2(0- 
The eonverted DAE is 

' 0 = 7i = --^1 - ^1 (0 - ^ 2(0 - 4 (0 + 4(0 + <^ 2(0 

_2 0 = /2 = -x'2 +X 4 + b2it) 

^ ■ { - 

0 = 4 = -XI -X2-Ciit)EC2it) 

, 0 = /4 = -.ri+ji;3+JC4 + C2(0- 




X] 

2^2 

2(^3 

X 4 

Ci 


Xl 

2C2 

2f3 

X 4 

/l 


0 * 

— 

— 

— 

1 

fl 

' -1 

0 

0 

0 

fl 


— 

1 

— 

0 * 

0 

fl 

0 

-1 

0 

1 

fl 


0 

0 * 

— 

— 

1 


-1 

-1 

0 

0 

h 


0 

— 

0 * 

0 . 

0 

h 

0 

0 

1 

1 



1 

1 

0 

0 

Val(E2) 

= 0 

det(j2) = 

1 



The SA succeeds on this converted DAE and gives struetural index V 5 = 2. Sinee m°o and are 
nonzero eonstants, and are always equivalent. 

7.1.2 ES method 

Eor in (7.1), the ES method eannot eonvert it to a strueturally regular DAE. We illustrate this 
argument with one partieular ehoiee of / G L in eaeh iteration of the ES method, and do not explore 
all possible eombinations of sueh ehoiees. To handle the limitation of the ES method, further 
development is required, whieh is left as future work. 
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Let u = (1,-1,1,-1)^. Then = 0. Using ( 6 .1-6.4) finds 

L= 11,2,3,41, 5 = |L| = 4, /={l,2,3,4|, and C = maxc, = 0. 

^ ^ iei 

Assume we pick 1 = 3. Using (6.9), we introduce yj for each = {1,2,4}: 


(7.2) 



(7.3) 


From (7.3), we construct g/s in (6.11). Since / = 3, we can exclude g 3 and ^3 in the converted 
DAE; see Remark 6 . 6 . 

By (6.10) and from (7.3), we write 


Xi=yi+X 3 , 4 = 3 ^ 2 -JC 3 , and X 4 =y 4 -X 3 . 


In (4.10), we 


substitute for in 

yi +^(^3 fi 

y2-X3 x '2 fi 

y 4 -X 3 X 4 / 2 ,/ 3,/4 

The converted DAE is 

0 = /i = -yi+bi{t) 

0 = /2= yA-y2 + b2{t) 

0 = / 3 = X2+y4+C\{t) 

0 = /4 =-jci+y4 + C 2 (t) (7.4) 

0 = gi = -yi+x'i -X3 
0 = g2 = -y2+x'2+X3 
0 = g 4 = -y 4 +X 4 +X 3 . 
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Xl 


X2 X 3 X 4 


E = 


/i 

fi 

fs 

U 

^1 

gi 

gA 


0 « 

0 

r 


0 - 

0 * - 
0 0 * 

0 0 


y\ 

0 * 


^2 c,- 


- - - - 0 * 


0 


0 

0 * 


Val(E) = 1 


/i 

fi 

h 

fA 

^1 
gi 
gA L 


Xl 

0 

0 

0 

-1 

1 

0 

0 


JC2 

0 

0 

1 

0 

0 

1 

0 


2f3 M 

0 0 


0 

0 

0 

-1 

1 

1 


0 

0 

0 

0 

0 

1 


3^1 

-1 

0 

0 

0 

-1 

0 

0 


det(J) = 0 


^2 

0 

-1 

0 

0 

0 

-1 

0 


3^4 

0 

0 

1 

1 

0 

0 

0 


If M = (1, —1,1, — 1,0,0,1)^, then Jw = 0. We use ( 6 .1-6.4) again to find 
L= {1,2,3,4,7}, 5 =|L|= 5 , 7 = { 3,4,5,6,7 }, and C = 

Sinee 7 = 3,4 G L and 

J 3 -C = J 4 -C = 0-1 = -1<0, 


maxc; = 1 . 
iel 


eondition ( 6 . 8 ) is not satisfied, and thus applying the ES method to (7.4) does not guarantee a striet 
deerease in the value of E. 
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7.2 Artificially modified pendulum ModPendB 

From Pend (3.8), we construct a problem ModPendB by performing a linear transformation on 
the state variables: 


X 


1 1 0 


Zl 

y 

= 

0 1 1 


Z2 

A 


1 0 1 


Z3 


The resulting DAE is 

0 = /i = (Z1+Z2)"+(Z1+Z2)(Z3+Z1) 

0 = /2 = (Z2+Z3)" + (Z2 + Z3)(Z3+Zl) -g (7.5) 

0 = /3 = (zi+zif + izi + zsf-L^- 



Zl 

Z2 

Z3 

C/ 


Zl 

Z2 

Z3 

/l 

2 

2 

0' 

0 

/i 

1 

1 

o' 

E®= /2 

0 

2 

2 

0 

j0= /2 

0 

1 

1 

/3 

_ 0 

0 

0 . 

2 

/3 

_ 2a 

2(a + /3) 

2/3. 

dj 

2 

2 

2 

Val(EO) 

= 2 

det(jO) = 0 



Here a=zi+Z2 and /3 = Z2 + Z3. 


7.2.1 ES method 

If u = (1, —1,1)^, then = 0 . We apply the ES method, and (6.1-6.4) give 

L={l,2,3}, 5 =|L|= 3 , /={1,2,3}, and C = rnaxc,-= C3 = 2. 


Since w is a constant vector, picking any ZeL={l,2,3} gives an equivalent converted DAE. We 
show below the conversion for case 1 = 1. 

Since L\ {/} = {2,3}, we introduce variables W2 and W3 corresponding to Z 2 and 23, respec¬ 
tively. Using (6.10) gives 


Z2=Z2 


id2~C) 


= W2 + 


U2 (rfi-C) _ 


Z3 


_ (d2-C) _ 


= W3-h 


Ml 

M3 (rfi-C) _ 
Ml 


= W2-Zl 


W 3 +Z 1 . 


(7.6) 
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To perform substitutions, we first write the derivatives z", Z 2 and z" in /i and /2 explieitly: 

0 = /l = Zi +Z 2 + (zi +Z2)(Z3 +Zl) 

0 = /2 = Z 2 +Z 3 + {Z2 + Z3){Z3 +Zl )-g 

Then we 

substitute for in 


„// j/ jf f f 


^2 
.y/ I 


I Jf Jf f 

^3+^1 ^3 h 

VV2—Zl Z2 /3 

1^3+21 Z3 /3 

Taking (7.6) into eonsideration, we find the resulting DAE (with gi and yi removed as Z = 1) 

0 = /i= w"+ W2(2 zi+W3) 

0 = /2 = (W2+W3)"+(w2 + W3)(2zi+W3) -g 

0 = /3= W 2 + (W2 + W3)^— 

0 = g2 = -Z2 + W 2 -Zl 

0 = g3 = -Z3+W3+ZI. 


(7.7) 


E = 



Zl 

Z 2 

Z3 

W2 

W 3 



Z\ 

Z 2 

Z3 

W2 

W 3 

7i 

■ 0 

— 

— 

2 * 

o' 

0 

/i 

2w2 

0 

0 

1 

0 ■ 

72 

0* 

— 

— 

2 

2 

0 


2m 

0 

0 

1 

1 

73 

— 

— 

— 

0 

0* 

2 J = 

I 3 

0 

0 

0 

2(w2+m) 

2m 

g2 

0 

0* 

— 

0 

— 

0 

gi 

-1 

-1 

0 

0 

0 

g3 

. 0 

— 

0* 

— 

0. 

0 

g3 

1 

0 

-1 

0 

0 . 

dj 

0 

0 

0 

2 

2 

Val(E) = 2 



det(J) = 

-4L2 



Here = W 2 + W 3 . We use equation /3 = 0 to obtain det(J): 
det(J) = — 4 ( 2 w 2 + 2 w 2 W 3 + W 3 ) = —4L^ 7 ^ 0. 


Henee SA sueeeeds on (7.7). Because = 1 is a nonzero constant, (7.7) and (7.5) are always 
equivalent. 
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7.2.2 LC method 

We show below how to apply the LC method to (7.5). The resulting DAE is relatively eomplicated, 
and its equivalence to the original problem requires two conditions to be satisfied. 

Let = (a,/3,—1/2)^. Then = 0. Using (5. 3-5. 6 ) gives 

7® ={1,2,3}, 0° = O, and L° = {l,2}. 

Since = a and w® = /3 are not identically nonzero, the converted DAE is equivalent to (7.5) only 
if 7 ^ 0 for the I we pick. 

Assume that M° = a = zi+Z 2 7^0- We pick I = 1 and replace /i by 

/1 = “i/i + ^ 2/2 + “ 3 / 3 ^ 

= (Zl +Z2)/l + (Z2 + Z3)/2-/372 

= {Zl +Z2){ZI+Z2)” + (Zl +Z2)^(Z3 + Zl) + (Z2 + Z 3 ) (Z2 + Z 3 )" + (Z2 +Z3)^(Z3 +Zl) 

-g{Z2+ Z3) -(Z1+Z2)(Z1+Z2)^^ - (7 + Z 2 )^ -(Z2+Z3)(Z2+ Z 3 )" “ (Z2 + 4)^ 

= [(Zl +Z2)^ + (Z2 +Z3)^] (Z3 +Zl) —g(z2+Z3) — (z} +z7^ “ (z2 + 4)^ 

= L^(Z 3 +Zl) —g(Z 2 + Z 3 ) — (Zi +z7^ ~ (Z2 + 4)^- 
The resulting DAE is 

(0 = fi =L2(z3+Zl)-g(z2+Z3)-(z'i+Z2)^-(z2 + Z3)^ 

7^7 <! 0 = /2 = (Z2+Z3)"+(Z2 + Z3)(Z3+Zl)-g 
i 0 = /3 = (Zl+Z2)^ + (Z2+Z3)^-T2. 



Zl 

Z2 

Z3 

C, 


Zl 

Zl 

Z3 

fl 

1 

1 

1 ■ 

I 

fl 

'-2a' 

-2(a + /3)' 

-2/3'" 

= /2 

0 

2 

2 

0 

ji= fl 

0 

1 

1 

/3 

. 0 

0 

0 . 

2 

h 

2a 

2 (a + /3) 

2/3 . 

dj 

2 

2 

2 

Val(Ei) 

= 3 


det(ji) =0 



We use a and /3 to denote zi +Z 2 and Z 2 + Z 3 , respectively. Let also 7 denote Z 3 +zi. Note 
that (a,/3,7) are actually {x,y,X) in (3.8) —this notation is for simplicity only, and we shall not 
substitute a,/ 3 , 7 forzi +Z 2 , Z 2 +Z 3 , andz 3 +zi, respectively. That is, we do not use the ES method 
here. 

Matrix is still identically singular. We let = (a, 2a/3' — 2/3a', a')^. Then = 0. 

Using (5. 3-5. 6 ) again gives 

7^ = { 1 , 2 , 3 }, 0^=0, and = {2}. 
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Suppose 

^ = (zi+Z2)(4 + 4)-(z2 + Z3)(zi+4) ^0. 

Then the eonverted DAE is equivalent to We replaee fz by 
72 = m}/(+m1/2 + M^/^' 

= af{ + 2{aj5' — a' ^)f 2 + cc' 

= a(LV-^/3'-2aV-2/3T)+2(aj8'-a'/3)(^+/3r-^) 

+ 2a'(a'2 + a^ + j8'2 + j8^) 

= a(LV-^/3')+2(a/3'-a'j8)(j8r-^)+2a'(a'2 + j8'2) 

= (zi Ezi) L^(z 3 +zi) -^(z'l +Z 2 ) 

+ 2 (zi+Z2 )(z2+4) - (4+4)(z 2 + Z3) (Z2+Z3)(z3+Zl)-,? 

+ 2(Zi+Z2) (4+4)^+ (4+4)^ • 

The resulting DAE is 

0 = /i = L^(z 3 +Zi) — g(z2+Z3) — (zi +Z2)^ “ (4 + 4)^ 

0 = fl — (zi +Z 2 ) +z4 “,?(zi +z4 

+2 (zi+Z 2 )(z 2 +Z 3 ) - (z'l+ z4(z 2+Z3) (Z2+Z3)(z3+Zl)-g 

+2(zi+z4 (4+4)^+(4+4)^ 

. 0 = /3 = (Zl+Z2)^ + (Z2+Z3)^-^^- 
Zl Z2 Z3 Ci 

7i r 1 1 1*1 0 

E 2 = 72 1 * 1 1 0 

/ 3 1. 0 0 * oj 1 

dj I I I Val(E^) = 2 

Jaeobian is eomplicated, and we do not show it here. Its determinant is 

det(J^) = -4L^(zi +Z 2 ) [(^1 +Z 2 )(z 2 + Z 3 ) - (z 2 +Z3)(zi +Z 2 )] 

= -4aL^(a/3'-j8a') 7 0, 
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since we already assume 

= a = zi+Z 2 7 ^ 0 and = al3'— ^a'0. 

Therefore SA sueeeeds and gives struetural index V 5 = 1. 

Now we eonsider the case a/3' — /3a' = 0. Sinee 
0 = /z' = 2aa' +2/3/3' and a^O, 
we have 

0 = a/3'-/3a' = a/3' + /3 • (/3/3')/a = p'{a^ + /a = p'L^/a. 

So /3' = a' = 0. Sinee = (a,2a/3' — 2/3a', a')^ = (a,0,0)^, the first row in is identieally 
zero and the LC method is not applieable here. 

Hence, the DAEs and are equivalent under the eonditions 

= a = zi+Z 2 7 ^ 0 and 112 =/3'= Z 2 +4 7 ^ 0. 

7.3 DAE (6.5) and LC method 

We show below that applying the LC method to (6.5) does not convert it into a strueturally non¬ 
singular DAE, beeause the condition for the LC method is not satisfied. Reeall (6.5) and its SA 
result. 

0 = /i=xi+e-"i-" 2 ^ 2 +gi (0 

0 = /2 = Xi +X2x'2 -hX2 +g2(0- 



Xl 

■^2 

Ci 


Xl 

JC2 

/l 

' 1* 

2 ' 

0 

/l 

-M 

-PX2 


0 

1* 

1 

h 

1 

X2 . 


dj 1 2 Val(E) = 2 det(J) = 0 

Here jj. = lfu = (1,/i)^, then = 0. Using (5. 3-5. 6 ) gives 

/={l, 2 }, 0 = 0 , and L = { 1 }. 

Let I = 1 and replaee /i by 

7i = uifi + M 2/2 = /i+M /2 

= X\ + [J. + gi{t) + ll{xi +X2 x'2 -\-x\-\- g2{t))' 

= Xl + P + gl{t) + II {x\ -I-X2X2 + {x'2)^ + 2X2X2+g'2{t)') 

= + Slit) + P {l + x[ -I-X2X2 + {X2)^ -|- 2X2^2 +^2(0) ■ 
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The resulting DAE is 

0 = f I = Xi+gi(t)+e -f-jc 2 x '2 E (-’‘^ 2 )^ +g 2 ( 0 } 

0 = f2=Xi+X2X2+X2+g2(t)- 



X\ 

21:2 

Ci 


Xl 

2 C2 


' 1* 

2 ' 

0 

/i 

-ail 

—a 11x2 


. 0 

1* 

J = 

1 

fl 

1 

X2 _ 

dj 

1 

2 

Val(E) = 2 


det(J) 

= 0 


Here a =x[ Ex 2 x '2 + (^ 2 )^ + 2-^2-^2 + -?2(0 ji = 

The eonversion step does not reduee the value of signature matrix nor produee a nonsingular 
Jaeobian, because (5.5) is not satisfied: 

c7(xi,m) = o{xi,ii) = l = l-0 = di-6. 






Chapter 8 

Conclusion and future work 


We identified two types of struetural analysis’s failure. For the first type, the system Jaeobian is 
strueturally singular, and the failure is likely due to hidden symbolie eaneellations. One way to 
handle this is to perform symbolie simplifieations before applying SA. 

We foeused on dealing with the seeond type, where SA fails in a less obvious way. In this ease, 
the Jaeobian is not strueturally singular but is still identieally singular. We proposed two symbolie- 
numerie methods for eonverting a DAE with sueh singular Jaeobian to an equivalent DAE on whieh 
SA sueeeeds with nonsingular Jaeobian, provided some eonditions are satisfied. Sueh eonditions 
ean be eheeked automatieally. These eonversion methods provably sueeeed and thus allow SA to 
handle more DAE types. Our methods enable SA to better reeognize the true strueture of a DAE, 
and thus SA is more likely to sueeeed and obtain eorreet struetural information. Moreover, our 
methods provide insights into reasons for SA’s failures, whieh were not well understood before. 

We summarize the two eonversion methods here. The EC method is more straightforward: it 
keeps the size of the system and replaees only one equation within a eonversion step. The ES 
method requires more eonditions to apply. It augments the system and ehanges several equations 
within a eonversion step, whieh generally takes more symbolie operations. The eommon goal of 
both methods is to reduee the value of the signature matrix; this value is also the number of degrees 
of freedom reported by the SA. We also need to ensure that the eonverted DAE is equivalent to 
the original one: on some time interval, a solution of the original DAE should be a solution to the 
eonverted DAE, and viee versa. Moreover, it is desirable to ehoose a eonversion (if possible) sueh 
that we do not need to monitor the equivalenee eondition (w/ 7 ^ 0 ) when we solve the eonverted 
problem. 

A praetieal question worth eonsidering is how to ehoose the appropriate eonversion method 
between the two for a given strueturally singular DAE. Eor many of the examples we have stud¬ 
ied, it is fairly eommon that eonditions of only one method are satisfied; see §7.1 and §7.3. Eor 
some other DAEs, provided eonditions of both methods are satisfied, applying one method usually 
requires fewer symbolie manipulations than applying the other; see ModPendB in §7.2. In the 
latter ease, we examine if a eonversion has an identieally nonzero ui, sueh that the eonverted DAE 
and the original one are always equivalent. If applying either method guarantees sueh equivalenee, 
we prefer the EC method beeause it ehanges only one equation and maintains the problem size. 
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Our next goal is to combine these conversion methods with block triangular forms (BTFs) of 
a DAE [28,29]. If the Jacobian is identically singular and the DAE has a non-trivial BTE, that is, 
it has two or more diagonal blocks, then we can locate the block that leads to the singularity, and 
perform a conversion step on this singular block. Using this approach, we have already made some 
progress and managed to convert several problems, including the Campbell-Griepentrog Robot 
Arm and the Ring Modulator, into SA success cases. However, more careful studies are required 
and details need to be worked out. 

Euture work also includes rigorous implementation of these conversion methods in MATLAB. 
We currently have a prototype code, which builds upon our structural analyzer package DAES A [23] 
and takes advantage of matlab’s symbolic toolbox. We can call daesa functions to display a 
converted DAE’s structure, perform quasilinearity analysis, and print out a solution scheme [24]. 
Our goal is to eventually build a complete tool for performing DAE conversions. 

Another interesting direction for research is combining the dummy derivative index reduction 
method [17] with the conversion techniques. McKenzie et al. [19] show that dummy derivative 
method and the E-method share some similarities. The former method converts a high-index DAE 
into a low-index one by introducing dummy algebraic variables and augmenting the system. Sup¬ 
pose we adopt this methodology. When a condition for applying a conversion method is violated, 
applying some dummy-derivative-like strategy may help make a conversion possible. 

Our conversion methods require a symbolic toolbox that is capable of performing nontrivial 
symbolic arithmetic operations, differentiation, and simplifications. However, we only focus on a 
linear combination of the equations (in the EC method) and a linear combination of derivatives of 
highest order (in the ES method). How to further exploit symbolic tools to develop other conversion 
methods should deserve some investigation. 

We conclude with a conjecture here. In all our experiments, when we transform a DAE with 
identically singular Jacobian to an equivalent solvable DAE with generically nonsingular Jacobian, 
the value of the signature matrix always decreases. As Pryce points out in [26], the solvability of 
a DAE lies within its inherent nature, not the way it is formulated nor the method that analyzes 
it. We conjecture that, if a reformulation of a structurally singular DAE results in an equivalent 
solvable DAE, then the value of the signature matrix always decreases. However, based on our 
current knowledge, it seems difficult to prove this conjecture. 


Appendix A 

Proofs for the ES method 


For readers’ eonvenience, we recall the notation from Sections 6.1 and 6.2: 


^={2 I “;^0} = {1:5}, ^=1^1, 

7 = { z I — Cj = Oij for some 7 G L }, C = rnaxc,. 

We also assume the conditions for applying the ES method are satisfied: 


o {xj,u) < 


dj-C-l if; gL 
dj — C otherwise, 


and 


(A.l) 


dj — C>0 for all j G L. 

Denote 

7 = {1 : zz} \ 7 = { z I dj — Ci> Oij for all j G L } , and 
L = {l:zz}\L=| j I zz^=0} = {5+l: n}. 

In the following, we assume we pick a column index Z G L in a conversion step using the ES 
method. We also assume zz/ 7 ^ 0 for all t in some I C M. 


A.l Preliminary results for the proof of Lemma 6.4 


Lemma A.l Let I ^ L. If (A.l) holds, then for an r E L\ {/}, 


Hr (di—C) 
xj,yr + —x] 

Ul 


< 


dj-C-\ z/;GL\{/} = {1,...,Z-1,Z + 1,...,.} 
dj — C z/; G LU {/} = { Z, 5 + 1,5 + 2,... ,zz }. 


(A.2) 
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Proof. According to our assumptions at the beginning of §5, no symbolie eaneellation oeeurs in a 
strueturally singular DAE. Henee the formal HOD and the true HOD of a variable in a funetion are 
the same. Using (4.7) gives 


, , Ur (d,-C) 

o {Xj,yr + —x) 

Ul 


Idl-C) 


(A.3) 


< max I (7 {xj, u ), o \ Xj,xi 

(a) Consider the ease j = I e L. Using (A.l) gives o{xi,u) < — C — 1, so 

RHS of (A.3) = max |(7 (jc/, m) , o j | 

= C7 (xi,xf'~^^^ =di-C. (A.4) 

(b) Consider the ease j 7 ^ /, that is, y G { 1,..., / — 1, / + 1,..., n }. Using (A. 1) again, we have 


RHS of (A.3) = max I <7 [xj, u ), o 


Xj,Xl 


{di-C) 


= <7 {Xj,u) < 


dj-C-l if; gL\{/} = { 5 } 

dj — C if ^ L = {/,5+1,... }. 


(A.5) 


Combining (A.4) and (A.5) results in (A.2) and eompletes this proof. 
Corollary A.2 For an i G /, 


□ 


f Ur {d,-c)\^c-ci)\ ( dj-Ci-\ (/■;GL\{Z} = { 1 ,. ..., 5 } 

a\xj,[yr+-x, ) 


otherwise. 


(A.6) 


Proof. Sinee C = max/g/c,, the order C — c,- > 0 for i G /. We have 
LHS of (A.6) = C7(.r;, 


j,yr + —xf^ ^^^+(C-q) 
Ul ) 


<(C-c,) + 


dj-C-\ if;GL\{Z} 


dj-C 


otherwise 
dj — Ci — 1 if;GL\{Z} 


using (A.2) 


dj-Ci 

= RHS of (A.6). 


otherwise 


□ 
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A.2 Proof of Lemma 6.4 

Proof. We write E from Figure 6 .1 into the following 2x3 bloek form 

E = 


^ 1,1 

^1,2 

^1,3 

Ml 

to 

7^2,2 

Ml 

to 


We aim to prove the relations between Oij and dj — Ci in eaeh bloek. 
Reeall (6.13) and (6.14): 

dj = dj, Ci = Ci if i,j=l:n 

dj = Ci = C if i,j — n +I : n + s. 

We prove for [Ei^ | Ei, 2 ], Ei, 3 , [E 2 ,i | ^ 2 , 2 ], and ^ 2,3 in order. 

1. Consider 


•^1 ■■■ a:/—1 Xi • • • Xs 'Tj-i-i ''' Xfi Ci 


"1,1 


" 1,2 


/l 


fn 


< 


< 


Cl 


dj di ■■■ di^i di ■■■ ds ds+i ■■■ dn 


We need to show that, for i 



= I : n, 

dj -Ci-\ 
dj - Ci 


ify GL = {1 :5}, 
if 7 G L = {5 + 1 : n}. 


We eonsider the cases 

(a) 7 G L and i G I, 

(b) 7 G L and i G I, and 

(c) 7 G L U L and i G 7. 

Let rGL\{/} = {l,— 1,Z+1,..., 5 }. 


(A.7) 
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(a) Consider J G E\{1}. By Corollary A.2, the HOD of Xf is < dj — c/ — 1 in every 


Ui 


(C-Ci) 




that replaces in an fi —here r G L\{ / }, which includes j. Therefore, 

a {xjJi) <dj-Ci-\ for j eL\{l}, i el. (A.8) 

Now consider j = I e L. We show below that dfjdx^f' = 0, which implies xf^' 
does not appear in /,, i G I. That is, 


(y {xiJi) <di-Ci-\ for i G I. 


(A.9) 


Using (A.l) gives 


o 



< o{xi,u) < di —C — \ . 


Also 


o 


Xl,yr + 


Ui ' J 


max |(7 (;c/, m) , o (^xi,x[ 


(di-C) 


= di-C. 


Since C — c,- > 0 for i G I, we apply Griewank’s lemma (5.1), with q = C — ci, to 


Ul 


(A. 10) 


from (6.10). Differentiating both sides of (A. 10) with respect to x 


{d,~c) 


grves 


Ur 

Ul 




^xidr-c+c-a) 




dx 


(di-C) 


dx 


(di—C+C—Cr 

I 


dx 


{di-Cj) 


(A.ll) 


Then 


dfi 

dxf~‘^‘^ 


^ fi 

rehl} 

h+ I J.-- 

reAlO 


/ . J iridr 


— {Ju)i = 0 

Ul 




by the chain rule 
by (A.ll) 


because Jw = 0. 
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(b) Now 7 G L = {5 + 1 : n} and i G I. None of the with j E L and i G / is replaeed 

By Corollary A.2, the HOD of xj is < dj — Ci in every 




—r'i\ ‘'d 


that replaees Henee 

<y (xj, fi) <dj- Ci for j ^L,iE I. 


(A.12) 


(e) Consider i G I. Sinee dj — Ci > Oij for every j G L, we do not replace any derivative. 
Hence, 


a{xj,fi)=o{xj,fi) < 


dj — Ci — 1 for / G /, 7 G L = {1 : 5 } 
dj — Ci for i G 7, 7 G L = {5 + 1 : n}. 


(A.13) 


Since dj = dj and Ci = q for i, j = \:n, combining (A. 8 , A.9, A.12, A.13) proves (A.7). 
2. For 


yi yi-\ yi yi+\ ■■■ ys Ci 


^ 1,3 = 


/i 


fn 


< 


< 


C\ 


di C ■■■ C C C ■■■ C 


we need to show that, for z = 1 : n, 

_ z if; eL\{Z}, 

Ct;>+7 — <7 [yj,fi) i -f ■ 1 

[ =-°o if7 = ^- 

Here, the {n + 7 )th column corresponds to yj. 

(a) Consider 7 G L\{1}. For all i G /, 

= 0 + (C-c,) = C-Ci. 
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Then for all i. 


( C- Ci Hi El and 




c<) is replaced by ("jj + 

J \ Ui 


— oo otherwise. 

To combine these two cases, we can write 

'^i.n+j ^C — Ci = dn+j — Ci for 7 G L\ { / } and all z = 1 : n. 

(b) Now consider j = 1. Since yi does not appear in any 

^ i,n+l Cr 

3. Consider 


- 2,1 


- 2,2 


fi)- 

— OQ 

for all z = 

1 : n. 







2Cl ■■■ 2C/_i 

Xl 

Xl+l ■■■ Xs 

■■■ 


Ci 

fn+l 

gl 

= 

= 


< 


C 


< 









< 




f n+l 

gl 


= 


— 00 . . . - 

— 00 

c 



< 


< . 









< 



fn+s 

gs 

- 

= 

= 


- 

c 


dj 

di ■■ ■ di-i 

di 

di+i ■■■ ds 

d-s+i 

d-n 



Recall I eL. Let row number z G {1 : 5 }. We consider the following cases. 

(a) j = l or j = i. That is, the entries in the Ith column in E 2,1 or those on the (main) diagonal 
of E 2 ,i. 

(b) j ^ I and z = 1. That is, the /th row in [E 2 ,i, E 2 , 2 ] without the Ith column. 

(cl) 7 = 1:5 and 7 , /, z are distinct. This case covers all the entries with ‘<’ in E 2 ,i. 

(c2) j = s + l : n and 7 , IJ are distinct. This case covers all the entries with ‘<’ in E 2 , 2 - 
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We shall prove 


— dj Cfi+i 




— dj Cfi-^i 1 

^ — dj Cn+i 


if 7 = ^ or j = i 

if j I and i = I 

if 7 = 1 : 5 , and 7 , /, i are distinct 

if 7 = 5 + 1 : n and 7 , 1, i are distinct. 

(A. 14) 


Recall 


0 = g/ = 


-yi+^ 


(di-C) 


Uj 


Ui 


.(di-C) 


■yi+x 


(di-C) 


for i G L\{ 1 } 
for i = 1. 


(a) Sincex 


and x^f' (if / = i then both are the same) occur in gi. 


o{xi,gi) = di-C = di- Cn+i and 
o{xugi) =di-C = di - Cn+i- 


(A.15) 


(b) Now 


cy = <7 [xj,yi-x\ 


Idi-C) 


= —00 <dj — C— \ = dj — Cn+l — 1 ■ 


(A.16) 


(cl, c2) Consider the last two cases together: 7 , /, and i are distinct. We have 
<y {xj,gi) = <y (^xj,yi-x\'^‘~^^ + <<y{xj,u). 

Using (A.l) gives 


(y {xj,gi) < (j {xj,u) < 


dj C 1 — d j CyiJ^i 1 


dj C 


— dj Cfi-\-i 


if 7 G L 
if 7 G I. 


(A. 17) 


Combining (A.15-A.17) gives (A.14). 
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4. For 


yi ■■■ yi^i yi yi+i ■■■ ys 


f n+l i?l 


£ 2,3 — f n+l Si 


f n+s Ss 


— 00 


0 


-00 


c 


c 


c 


di C ■■■ C C C ■■■ C 


we consider z, j = I : s. Then 


(yn+i,n+j = (y {yj ,gi)=o[yj, yi - ^ 


(di~C) ^ Ui (d,~C) 


= (yjAi) 

0 C C dfi+j Cn+; 


if i = j 
otherwise. □ 



Appendix B 

More examples 


We review how to remedy several strueturally singular DAEs from the literature. We diseuss the 
index-5 Campbell-Griepentrog Robot Arm DAE in §B. 1 , the transistor amplifier DAE in §B.2, and 
the ring modulator problem in §B.3. We also show in §B.4 how to “fix” the index overestimation 
problem on ReiBig’s family of linear DAEs, for whieh SA produees a nonsingular system Jaeobian 
but overestimates the index. After we apply a teehnique similar to the linear eombination method, 
SA reports the eorreet index V 5 = 1 on this family of DAEs. 


B.l Robot Arm 


We slightly simplify the two-link robot arm problem in [5] by writing the derivatives of xi,X 2 , and 
X 3 implieitly in the equations: 


0 = /i 

0 = /2 

0 = /3 

0 = /4 
0 = /5 


ff 

- 


ff 

^2 


2c (X3 ) {x\ -|- X3 ) ^ -|- x'^d (X3 ) 

-b (2X3 -X2){a{x3) + 2b{x3)) +a{x 3 ){ux-U 2 ) 

- -2c (X3) {x\ + X3)^ - x'ld (X3) 

-|- (2x3 ~ — 2b{x3)) — a{x3)ui + («(X3) -|- 1)m2 

x" - -2c(x3) (x'l -fXg)^ -xf <i(x3) -b (2X3 -X 2 ) (<3(x3) - 9 b{X3)) 

-2xfc(x3) -<i(x3)(xi -fXg)^- (a(x3) +b{x 3 )){ui - U2) 

eosxi -|-eos(xi -I-X3) —pi{t) 
sinxi -|-sin(xi -I-X3) —p2{t), 


(B.l) 
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where 


piW = 

cos(l 

— e^) +eos 

(1- 

-t) 

P2{t) = 

sin(l —e 

+sin(l -t) 


a{s 

) = 

2/(2- 

2 

-COS 




b{s) = 

coss/ (2 

-eos^i) 


c{s 

) = 

sini/(2 — eos^i) 



d{s) = 

sin^eosi 

/(2 —cos^i). 



X\ 

-^2 

^3 

Ml 

M2 

Ci 


Xi X2 

X3 Ml 

M2 

fi 

■ 2 

0 

0 

0* 

o' 

0 

fi 

' 1 

-«3 

«3 

fi 

s 

2* 

0 

0 

0 

0 

h 

1 

<33 

-1 -(23 

h 

s 

0 

2 

0 

0* 

0 

J= /3 


1 <23+^73 

-(23 - b-i 

/4 

0 


0* 



2 

/4 

df4 

c)xi 

dh 

dx^ 


/s 

.0* 


0 



2 

h 

dh 

-dx\ 

dh 

dx^ 


dj 

2 

2 

2 

0 

0 

Val(E) 

= 2 


det(J) = 0 



Here 

df^jdxi — — sin;ci — sin(;ci +X3) = a{x2,) = Ij (2 — eos^X3) 

df^jdx^ — — sin(xi +;C3) = b{x2,) = eosji;3/(2 — eos^X3) 

df^/dxi = cosxi +eos(xi +^3) 
dfs/dxs = cos(a:i +X3). 

SA reports struetural index 3, while the d-index is 5. 


B.1.1 ES method 

Pryee [26] fixes this failure by introdueing a new variable w and substituting it for mi — M 2 in /i 
and / 3 . These two equations beeome /j and We append g = w — (mi — M 2 ) that preseribes these 
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substitutions and obtain the converted DAE: 

0 = fi =Xi — 2 c(x 3) {x'l -\-x!^)^ + x!id{x3) 

+ {2X3 -X2){a{x3)+2b{x3)) +a{x3)w 
0 = /2 =X2 - -2c{x3) (/i + 4 )^ -xfd{x3) 

+ {2x3 —X 2 ) (1 — 2>a{x3) — 2b{x3)'j —a{x3)ui + {a{x3) + 1 ) m 2 
0 = /3 =X 3 - -2c{x3){x\ +x' 3 f -x'^d{x 3 ) + { 2 x 3 -X 2 ) (^(jcs) -9b{x3)) 
- 2xfc{x3)-d{x3){x\+X3)^ - {a{x3)+b{x3))w 

0 = /4 =cosa:i +cos(a:i +X3) — pi{t) 

0 = /s =sina:i +sin(a:i + 2 ^ 3 ) -P 2 {t) 

0 = g=w-{ui-U 2 ). 



Xl 

-^2 

-^3 

Ml 

U2 

w 

Ci 

/l 

' 2 

0* 

1 



o' 

2 

/2 

1 

2 

1 

0‘ 

0 


0 

/s 

1 

0 

2 



0* 

2 

/4 

0 


0* 




4 

/s 

0* 


0 




4 

g 

_ 



0 

0* 

0. 

0 

dj 

4 

2 

4 

0 

0 

2 

Val(E) = 0 



2Cl 

X2 

X3 Ml 

M2 W 

7i 

' 1 

«3 + 2^3 


-<33 

/2 


1 

<33 

-<33 - 1 

73 


<33 - 9b3 

1 

< 33+^3 

/4 

dfA 


df4 


()x\ 


dx^ 


/s 

ill 


ill 


dxi 


dx^ 





-1 

1 


det(J) = —2smx3{al — 3a3b3+b2) 


(B.3) 
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J is not identically singular; refer to (B.2) for the entries in it. SA reports the correct index 5, and 
succeeds if det(J) ^ 0. 


B.1.2 LC method 

We replace fs by 


/s — /i + 


as 


-fs 


= JC, 


+ 


as+bs 

2 c(xs)(x'i +X3)^+xf J(X3) + (2X3 -Jf 2)(«3 +2^3) + < 33(^1 -^2) 
as 


as // 

■•*^3 ~ ■ 


as+bs as + bs 


2c(x3)(xi +X2y-Xid{xs) + (2x3 -m) {as- 9 bs) 


ff 

= Xl - 


- 2 xf c(x3) -d{xs) (x'l +X 3 )^ -(a 3 + h3)(Mi - M 2 ) 
2 c(x 3 )(Xi +X3)^+xf^i?(x3) + ( 2 X 3 “-^2) (<33 + 2^3) 


+ -^4 

as+bs 


as 


- 2 c{xs){x\ +x'2,y -xfd{xs) + ( 2 X 3 -X2) {as- 9 bs) 


as+bs L 
2 xf c(X3 )-d{xs) (x 1 + X3)' 


(the underlined terms cancel out). 



XI 

4:2 

4:3 

Ml 

M2 

Ci 


XI 

X2 

xs 

Ml M2 

/l 

r@ 

0 

0 

0* 

o' 

0 

fl 




1- 

1 

fl 

□ 

2 * 

0 

0 

0* 

0 

fl 


1 


1 

7 

(T^ 

Is 

2 

0* 

2 



2 J = 

fs 

1 

fh 

dx 2 

03+^3 


/4 

0* 


0 



4 

h 

dfA 

6 x 1 


dfA 

3x2 


/s 

. 0 


0* 



4 

fs 

-dx\ 


dxj, 

- 

dj 

4 

2 

4 

0 

0 

Val(E) = 0 







Here 


dfs 


= as+ 2 bs + 


as 


-{as- 9 bs), 


(9x2 «3 + ^3 

det(J) = —2 sinx 3 (a 3 — 3 ( 33 ^ 3 +h 3 )a 3 /(a 3 +h 3 ). 
Refer to (B.2) for the other entries in J. Since 
as 


as+bs 2 + COSX3 


7 ^ 0 for all X 3 G 
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the converted DAE is always equivalent to (B.l). J is not identically singular. SA reports index 5, 
and succeeds if det(J) 7 ^ 0. 

B.2 Transistor amplifier 

Below is a transistor amplifier problem originated from electrical circuit analysis [18]. It is classi¬ 
fied in [18] as a stiff index-1 DAE consisting of 8 equations. 


0 = /5 


0 = /7 


0 = /6 


0 = /4 


0 = /2 


0 = /l 


0 = /3 



(B.4) 


0 = /8 


Cs{xq—x^) E 


where 


g{y)= I5{exp{y/Uf)-1) 


Ueit) = 0.1 sin( 200 ;rr) 


Ub = 6.0 
Up = 0.026 
a = 0.99 
/3 = 10 ^ 


Rq = 1000 

Rk = 9000 for/t= 1,...,9 

Ck = kx for/t= 1,...,5 





APPENDIX B. MORE EXAMPEES 


70 


Xi X 2 X 2 , X4 X 5 X^ X-j Xs c, 

/i r 1 1* 10 

/2 1* 1 0 0 

/3 0 !• 0 

/4 0 0 1 r 0 

/5 !• 1 0 0 
/6 0 r 0 

/7 0 0 1 r 0 

h [ r ij 0 

d; 1 1 1 1 1 1 1 1 Val(E) = 8 


J = 


Xl X 2 XT, X4 X5 

h r Cl -Cl 

fi -Cl Cl 

/3 C2 

/4 C 3 -C 3 

/s -C 3 C 3 

/6 
fl 

h . 


X6 


C4 


C5 

-C 5 



det(J) = 0 


SA reports index 1, but produces an identically singular J. Observing its structure, we 

replace by 
fl fl=fl+f 2 

f 4 74 =/ 4+/5 

fl fi = fiEh- 
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The new equations in the converted DAE are 

Xi-Ue(t) t/. / 1 1 \ , , 

0 = /a = ~^3) - ^ +^5 + -(a-l)^(;c5-;c6) 

- xj-Ub . x^ 

^ = fi = s — + (^si^5-xe) + -^■ 

Kg Kg 


-^1 

-^2 

X 3 X 4 

-^5 

;C6 

-3:7 

0* 

0 

0 




1 

r 

0 





0 

1* 





0 

• 

0 

0 

0 

0 




1 

r 

0 





0 

r 





0 

0 

0* 






1 

1 

1 

1 1 

1 

1 

1 

x\ 

-^2 

-^3 

;C4 

-^5 

X6 

Ro' 

7x2 

7x3 




-Cl 

Cl 

C2 





7/4 

dfl 

7/4 

7x3 

^4' 

7/4 

7x5 

7/4 

7x6 




-C3 

C3 



dfl dfl 
7x5 7x6 


R-^' J!,-‘ 


det(j) = C1C2C3C4C5 k-' + fc' + (;(8 ‘ 


-C5 C5 

■H) (''d 
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Note that only two partial derivatives shown in J contribute to det(J): 


df4/dx5=R5^+R^^ 

dfi/dx 2 =Ri^+R 2 ^- 


SA still reports index 1. Since J is not identically singular, SA succeeds if det(J) ^ 0. 


B.3 Ring modulator 


Following is a ring modulator problem originated from electrical circuit analysis [18]. By setting 
= 0 in the original problem formulation, we obtain an index-2 DAE consisting of 11 differential 
and 4 algebraic equations: 


0 = /i =-y[+C ^(y8-0.5yio + 0.5yii-Fyi4-i? Vi) 

0 = /2 = -y2-FC^^(y9-0.5yii+0.5yi2+yi5-^^V2) 

0 = /3 = yio —‘I(Udi)+<i(Ud4) 

0 = /4 = —yn + (i(Ud2) — (i(Ud3) 

0 = /5 = yi 2 + g(UDi)-gi^Ds) 

0 = /6 =—yi3 — g(UD2)+g(UD4) 

0 = /7 = -yj+Cp^-Rp^yj+giUm) +g(UD2) - g(UD3)-q(UD4)) 

0 = /8 =-ys + -Lj;^yi (B.5) 

o = /9 =-y9 + -L/;^y2 
0 = /io = ->'io + ^ 72 ^( 0 - 5 ji ->^3 -%Jlo) 

0 = /ii = ->'ii+f^73^(-0.5yi+y4-^g3>'ii) 

0 = /l 2 = ->'12 + ^72^ (0.5y2 -J5 -Rg 2 y\ 2 ) 

0 = /l 3 = ->'l 3 +^ 73 ^(- 0 - 5 >' 2 +j 6 --^g 3 >' 13 ) 

0 = /l4 = —y'\4-^Rsl (~^l + Uin\(t) — (i?,- +i?gi)yi4) 

0 = /i5 = -y'i5+Rsi i~y2 - {Rc+Rgi)yi5), 


where 


Udi= y3—y5—yi — Uin2it) 
Ud 2 = —y4 Eye —yi — ^i«2(0 
Ud 3= y4+y5+yi EUin2{t) 
Ud4 = —J3 —ye+yi EUin2{t) 


q{U) = y{e^^-\) 
Uin\{t) = 0.5 sin(2000;rt) 
Ui„ 2 {t) = 2sin(20000;rt) 
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The parameters are 

C= 1.6x 10^^ 
Cp = io-« 

R = 25x 10 ^ 
Rp = 50 
Lh = 4.45 

Lsi=2x 

L,2 = 5 X 10-4 
L,3 = 5 X 10 - 4 . 


Rgi = 36.3 
Rg2 = 17.3 
Rg3 = 17.3 
Ri = 5x 10 
Rc = 6x 10^ 

7 = 40.67286402 X 10-‘^ 
5 = 17.7493332 


SA reports index 1 and produees the following E with Val(E) = 11. 


yi yi y 3 y 4 ys ye yi y?, y 9 3^10 jn yn yu 3^15 


E = 


fi 

fi 

h 

/4 

h 

fe 

fl 

h 

/9 

/lO 

fll 

fn 

fn 

fu 

/l5 


r 


0 

0 


r 


0 

0 


0 0 0 * 0 

0 0 * 0 0 

0 * 0 0 0 

00 * 00 

0 0 0 0 r 


r 


r 


0 0 


r 


r 


0 0 


r 


r 


r 


dj 


1100001111 1 1 1 1 


We do not present J here. The entries that eontribute to its determinant are positions (/,z) for 
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i = 1,2,7 : 15, where dfi/dy\ — —1, and the submatrix 





3^4 

ys 

ye 


h 

-5l - 

54 

■^1 

-54 

J(3:6,3 

:6)= /4 


-52-' 

S3 -S3 

S2 


fs 

^1 

-■^3 

-Si - S3 



fe 

. -■^4 

^2 


-52-5 

Since this submatrix of J is identically singular, so is J. 

To remedy this DAE, we replace one equation from /s,/a,/s,/6 

/ — /s “/4 +/5 “/6 — JlO +yil +Jl2 +yi3- 


Consider the case /s. SA produces E with Val(E) = 10: 



yi y2 

y3 y4 

ys ye 

yv ys y9 

yio yii 

fi 

1* 



0 

0 0 

fi 

1* 



0 


h 





0 0 

U 


0* 

0 0 

0 

0 

/5 


0 0 

0* 

0 


/6 


0* 0 

0 

0 


fl 


0 0 

0 0 

1* 


Ml 

II 

0 



1* 


/9 

0 



1* 


/lO 

0 

0 



1* 

fll 

0 

0 



1* 

fn 

0 


0 



fn 

0 


0* 



fu 

0 





fl5 

0 





dj 

1 1 

0 0 

0 0 

1 1 1 

1 1 


, where 5,-= . 


0 

0 * 


r 


r 
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The consistent initial values for solving this converted DAE are 

yi = 0 for z = 1 : 15, 

y'. = 0 for z = 1,2,7 : 15, 

which satisfy ft for all i and /g. At this consistent point, det(J) = —1.2040 x 10^Hence SA 
succeeds and the DAE is index-2. 


B.4 An index-overestimated DAE 


ReiBig et al. [31] construct a family of linear constant coefficient DAEs 
Ax{t) +Bx{t) —q{t) = 0 


(B. 6 ) 


for which SA finds an arbitrarily high structural index V 5 > 1, though the d-index is 1. Here B is 
an zz X zz identity matrix, where n = 2k+\ and > 1; A has the form 


A = 


/Oil 
1 1 0 
0 1 1 
1 1 

0 


\ 


0 


V 


••• 0 
0 1 1 
1 1 
0 / 


That is, for i=l\k, equations f 2 i-\ and f 2 i have the common expression -^21 +-^ 2 i+i- 
Pryce [27] applies the E-method on (B. 6 ) with zz = 5 and k = 2\ 

0 = /i = X 2 +Xj,+Xi - qi (t) 

0 = /2 = X 2 +X 2 +X 2 - q2{t) 

0 = /3 = x'^+x'^+x-i - q 2 ,{t) 

0 = /4 = X 4 -hX5-hX4 — (2'4(t) 

0 = / 5 = x5-q5{t)- 


(B.7) 
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Xl 


^(^3 

X 4 

^5 

Ci 


Xl 

2 f 2 

^(^3 

X 4 

^5 

fl 

’o* 

1 

1 



0 

/i 

■ 1 

1 

1 



h 


r 

1 



0 

fl 


1 

1 



h 



0 * 

1 

1 

' J = 

h 



1 

1 

1 

/4 




1 * 

1 

1 

/4 




1 

1 

/s 





0 *. 

2 

/s 





1 _ 

dj 

0 

1 

1 

2 

2 

Val(E) = 2 



det(J) 

= 1 



The method sueeeeds with det(J) = 1 but reports V 5 = 3 different from = 1; this oeeurs in 
Pantelides’s SA as well. We illustrate below how to fix this index overestimation problem for 
(B.7). 

Observing the strueture of A, we 


replaee by 

/i 7 i=/i-/2 

/3 73 =/3-/4- 

The eonverted DAE is 


0 = /i =v:i-X2-<?i(t)+<?2(0 
0 = /2 = 4 + X 3 +X2 - q2{t) 

0 = f^=X'i-x^-q2{t)+qA{t) (B. 8 ) 

0 = /4 = X4 + V5 +X 4 — q 4 (t) 

0 = /5 =X5-q5(t). 


E = 



Xl 

2 C 2 


X 4 


Ci 


Xl 

X 2 X 3 

X 4 X 5 

7i 

’o* 

0 




1 

fl 

' 1 

-1 


fl 


1 * 

1 



0 

fl 


1 1 


fl 



0 * 

0 


1 

J= 73 


1 

-1 

/4 




1 * 

1 

0 

/4 



1 1 

fs 





0 *. 

1 

fs 



1_ 


Val(E) = 2 


det(J) = 1 


1 


1 


1 


1 


1 
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Since no dj = 0, SA reports index V 5 = max, c, = 1. Note that here we do not choose the canonical 
offsets 

c= ( 0 , 0 , 1 , 0 , 1 ) and d = ( 0 , 1 , 1 , 1 , 1 ), 
which still give an overestimated structural index Vs = 2 as =0 and C 3 = C 5 = 1 . 

Consider for general case k>\. The DAE is 

0 = fii -1 = X 2 i + 4/+ 1 + 1 - d 2 i~ 1 (0 i=\-.k 

0 = fii = 4; +4i+i+ Mi - quit) i=\:k (B.9) 

0 =f2k +1 = X2k +1 — q2k +1 (0 • 


E = 


Xi X2 V3 X4 X5 ■■■ X 2 k X2k~^l Ci 


/i 0* 1 1 

/2 1 * 1 

/3 0 * 

U 

f2k~l 

f2k 

f2k+l . 

dj 0 1 1 


1 1 

r 1 

0 * 

••• 1 

r 

2 2 • • • k-i 


0 

0 

1 

1 

1 k- 1 

1 k- 1 

0* J k 

k Val(E) = k 
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J = 


Xi X2 XT, X4 X5 ■■■ X2k X. 2 k+\ 

/i r 1 1 1 

fi 1 1 

h 111 

/4 1 1 


flk-l 

hk 

flk+l 


1 1 1 
1 1 
1 

det(J) = 1 


SA succeeds and reports structural index k+\. 

To remedy this index overestimation, we repeat the same strategy used above: for i = 
replace f 2 i~i by f 2 i~\ = fn-i — fa- The converted DAE is 

0 =f2i- 1 = 1 - X 2 i - i{t) A quit) i=\\k 

0 = f2i = X 2 i + X 2 i+ 1 + X 2 i -q2i{t) i=\:k 

0=/2/t+l =X2k+l — <?2/t+l(0- 


E = 


X\ X 2 XT, X4 X5 


fi 


0 * 


fl 

73 

h 


0 

r 1 
0 * 0 
r 1 

0 * 


flk-l 

flk 

flk+l 


rf; 1 1 1 1 1 


X2k X2k+\ c, 

1 

0 

1 

0 

0 1 
r 1 0 

0* J 1 

1 1 Val(E) = k 


1 : A:, we 


(B.IO) 
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J = 


Xi X2 X2 X4 X5 


fl \ 1 
/2 
/s 

/4 


-1 

1 1 
1 -1 
1 1 


flk-l 

flk 

flk+l 


Mk Mk+\ 


-1 

1 1 
1 


det(J) = 1 


Now SA reports the eorreet V 5 = 1 on the eonverted DAE (B.IO). Again, we use non-eanonieal 
offsets in E, while = ci = 0 in the canonieal ease. 
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